Taller 1 - Datos funcionales¶
Librerías¶
import os
import warnings
import numpy as np
import scipy.io as sio
import pandas as pd
import matplotlib.pyplot as plt
import skfda
from skfda.representation.basis import BSplineBasis
from skfda.preprocessing.smoothing import BasisSmoother
from skfda.exploratory.stats import mean, trim_mean, cov, depth_based_median
from skfda.exploratory.depth import ModifiedBandDepth, BandDepth
from tqdm import tqdm
from skfda.exploratory.visualization import Boxplot
# Suppress warnings
warnings.filterwarnings("ignore")
Cargar los datos¶
# Load MATLAB data
archivo_mat = "../Datos/data.mat"
datos = sio.loadmat(archivo_mat)
# Reshape data
# Dimensions: 268 tests x 571 measurement points x 7 wavelengths
data = datos["X"].reshape(268, 571, 7, order='F')
# Wavelength labels
wavelengths = {
0: 'Y1=230', 1: 'Y2=240', 2: 'Y3=255',
3: 'Y4=290', 4: 'Y5=305', 5: 'Y6=325', 6: 'Y7=340'
}
# Measurement points grid
dominio = np.arange(275, 560.5, 0.5)
# longitud de onda a trabajar
wavelength_idx = 6
Visualization of Raw Data¶
def plot_wavelength_data(organizado, n_samples=100):
"""
Plot intensity data for first n_samples across 7 wavelengths
"""
colors = plt.cm.tab10(np.linspace(0, 1, n_samples))
fig, axes = plt.subplots(4, 2, figsize=(20, 25))
fig.suptitle('7 Wavelengths - First 100 Samples', fontsize=16)
axes = axes.ravel()
x = np.arange(571)
for wavelength in range(7):
for sample in range(n_samples):
data = organizado[sample, :, wavelength]
axes[wavelength].plot(x, data, color=colors[sample],
linewidth=1.5, alpha=0.8)
axes[wavelength].set_title(f'Wavelength {wavelength+1}', fontsize=12)
axes[wavelength].set_xlabel('Measurement Points')
axes[wavelength].set_ylabel('Intensity')
axes[wavelength].grid(True, alpha=0.3)
axes[-1].set_visible(False)
plt.tight_layout(rect=[0, 0.03, 0.95, 0.95])
plt.show()
plot_wavelength_data(data)
1.1 Suavizado spline¶
def smooth_functional_data(
data,
wavelengths,
dominio,
base_path="../Resultados/punto_1",
n_basis=65,
best_param=1e-3,
rango_curvas=[None, None],
mostrar_imagen_cada=None
):
"""
Perform functional data smoothing for specified range of samples and wavelengths
Parameters:
-----------
data : numpy.ndarray
Input data array with dimensions (n_samples, n_points, n_wavelengths)
wavelengths : dict
Dictionary of wavelengths
dominio : array-like
Domain points for the data
base_path : str, optional
Base directory to save results
n_basis : int, optional
Number of B-Spline basis functions
best_param : float, optional
Smoothing parameter (lambda)
rango_curvas : list, optional
Range of curves to process [start, end].
If None, processes all curves.
Indices are 0-based.
mostrar_imagen_cada : int or None, optional
Show an image every 'n' processed samples.
If None, no images are shown during processing.
"""
# Crear la estructura de carpetas base
if not os.path.exists(base_path):
os.makedirs(base_path)
# Crear subcarpetas para cada longitud de onda
for wave_idx, wave_name in wavelengths.items():
wave_path = os.path.join(base_path, wave_name)
if not os.path.exists(wave_path):
os.makedirs(wave_path)
# Determinar el rango de curvas a procesar
inicio = rango_curvas[0] if rango_curvas[0] is not None else 0
fin = rango_curvas[1] if rango_curvas[1] is not None else 5
# Crear una barra de progreso
total_iterations = fin - inicio
with tqdm(total=total_iterations, desc="Procesando muestras") as pbar:
# Contador de muestras procesadas
muestras_procesadas = 0
# Para cada longitud de onda
for wave_idx, wave_name in wavelengths.items():
# Si ya hemos procesado todas las muestras requeridas, salir del bucle
if muestras_procesadas >= fin:
break
# Procesar muestras para esta longitud de onda
for muestra in range(data.shape[0]):
# Si ya hemos procesado todas las muestras requeridas, salir del bucle
if muestras_procesadas >= fin:
break
# Obtener los datos para esta muestra y longitud de onda
X = data[muestra, :, wave_idx].flatten()
# Crear una base B-Spline
basis = BSplineBasis(domain_range=(275, 560), n_basis=n_basis)
# Crear functional data grid
fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio)
# Suavizado penalizado
smoother = BasisSmoother(basis, smoothing_parameter=best_param)
fd_smoothed = smoother.fit_transform(fd)
# Decidir si mostrar la imagen
if mostrar_imagen_cada is not None and muestras_procesadas % mostrar_imagen_cada == 0:
# Crear la figura
plt.figure(figsize=(12, 7))
# Graficar los datos originales como puntos
plt.scatter(dominio, X, c='blue', s=30, alpha=0.5, label='Datos originales')
# Graficar la curva ajustada
plt.plot(dominio, fd_smoothed(dominio)[0], 'r-',
label=f'Ajuste B-spline',
linewidth=2)
# Personalización del gráfico
plt.title(f'Longitud de onda {wave_name} nm - Muestra global {muestra}\n'
f'Bases B-spline: {n_basis}', fontsize=12)
plt.xlabel('Longitud de onda (nm)', fontsize=11)
plt.ylabel('Intensidad', fontsize=11)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
# Ajustar los márgenes
plt.tight_layout()
# Mostrar la figura
plt.show()
# Crear la figura para guardar
plt.figure(figsize=(12, 7))
# Graficar los datos originales como puntos
plt.scatter(dominio, X, c='blue', s=30, alpha=0.5, label='Datos originales')
# Graficar la curva ajustada
plt.plot(dominio, fd_smoothed(dominio)[0], 'r-',
label=f'Ajuste B-spline',
linewidth=2)
# Personalización del gráfico
plt.title(f'Longitud de onda {wave_name} nm - Muestra global {muestra}\n'
f'Bases B-spline: {n_basis}', fontsize=12)
plt.xlabel('Longitud de onda (nm)', fontsize=11)
plt.ylabel('Intensidad', fontsize=11)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
# Ajustar los márgenes
plt.tight_layout()
# Guardar la figura
output_path = os.path.join(
base_path,
wave_name,
f'onda_{wave_idx}_muestra_global_{muestra}.png'
)
plt.savefig(output_path, dpi=300, bbox_inches='tight')
plt.close()
# Incrementar el contador de muestras procesadas
muestras_procesadas += 1
# Actualizar la barra de progreso
pbar.update(1)
print("Proceso completado. Se han generado todas las imágenes con el ajuste penalizado.")
total_curvas = 7*268
print(f'Hay un total de {total_curvas} curvas a procesar')
smooth_functional_data(data, wavelengths, dominio, rango_curvas=[0, 5], mostrar_imagen_cada=1)
Hay un total de 1876 curvas a procesar
Procesando muestras: 0%| | 0/5 [00:00<?, ?it/s]
Procesando muestras: 20%|██ | 1/5 [00:01<00:06, 1.68s/it]
Procesando muestras: 40%|████ | 2/5 [00:03<00:05, 1.70s/it]
Procesando muestras: 60%|██████ | 3/5 [00:05<00:03, 1.72s/it]
Procesando muestras: 80%|████████ | 4/5 [00:07<00:01, 1.85s/it]
Procesando muestras: 100%|██████████| 5/5 [00:08<00:00, 1.77s/it]
Proceso completado. Se han generado todas las imágenes con el ajuste penalizado.
1.2¶
Estadistica descriptiva funcional¶
def plot_covariance_3d(cov_matrix):
"""
Visualiza una matriz de covarianza en un gráfico 3D de superficie.
Parámetros:
cov_matrix (numpy.ndarray): Matriz de covarianza a visualizar
"""
# Crear coordenadas para la malla
x = np.arange(cov_matrix.shape[0])
y = np.arange(cov_matrix.shape[1])
X, Y = np.meshgrid(x, y)
# Crear figura y eje 3D
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
# Graficar la superficie de la matriz de covarianza
surf = ax.plot_surface(X, Y, cov_matrix,
cmap='viridis', # Puedes cambiar el colormap
edgecolor='none',
alpha=0.8)
# Personalizar el gráfico
ax.set_title('Visualización 3D de Matriz de Covarianza', fontsize=15)
ax.set_xlabel('Dimensión 1', fontsize=12)
ax.set_ylabel('Dimensión 2', fontsize=12)
ax.set_zlabel('Valor de Covarianza', fontsize=12)
# Añadir una barra de color
fig.colorbar(surf, shrink=0.6, aspect=10)
plt.tight_layout()
plt.show()
def functional_stats_analysis(data, wavelength_idx):
"""
Compute and visualize functional statistics
"""
X = data[:, :, wavelength_idx]
fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio)
# Mean computation.
mean_func = mean(fd)
mean_trim_func = trim_mean(fd, 0.1)
# Plotting means
plt.figure(figsize=(12, 6))
for sample in range(268):
plt.plot(dominio, X[sample], color="blue", linewidth=1.5, alpha=0.3)
plt.plot(dominio, np.mean(X, axis=0), color='red', linewidth=3, label='Functional Mean')
plt.plot(dominio, mean_trim_func(dominio)[0], color='green', linewidth=3, label='Trimmed Mean (10%)')
plt.title(f'Individual Curves and Mean Functions - Wavelength {wavelengths[wavelength_idx]}')
plt.xlabel('Wavelength (nm)')
plt.ylabel('Intensity')
plt.legend()
plt.grid(True)
plt.show()
# Variance
var_func = skfda.exploratory.stats.var(fd)
var_func.dataset_name = f"Varianza Funcional"
var_func.plot(color='purple', linewidth=2, legend=True)
# Depth-based median
median_func = depth_based_median(fd)
median_func.dataset_name = f"Mediana basada en profundidad MBD"
median_func.plot(color='orange', linewidth=2, legend=True)
# covariance
covariance = cov(fd, correction = 0)
cov_matrix = covariance(fd.grid_points[0], fd.grid_points[0])
plot_covariance_3d(cov_matrix)
return fd
fd_analysis = functional_stats_analysis(data, wavelength_idx)
Profundidad y class_outliners¶
Analisis de profundidad funcional por BD y MBD¶
Definición de Band Depth (BD)¶
La Band Depth (BD) mide la "centralidad" de una función $x$ en relación con un conjunto de curvas $x_1, x_2, \ldots, x_n$. Se basa en calcular la proporción de bandas, formadas por subconjuntos de $j$ curvas, que contienen completamente a la función $x$.
Bandas Formadas por Curvas: Se tiene que el gráfico de la función $x$ es el subconjunto del plano definido por $G(x) = \{(t, x(t)) : t \in I\}.$ Para $j \geq 2$, una banda $B(x_{i_1}, x_{i_2}, \ldots, x_{i_j})$ delimitada por $j$ curvas $x_{i_1}, \ldots, x_{i_j}$ se define como:
$$ B(x_{i_1}, x_{i_2}, \ldots, x_{i_j}) = \{ (t, y) : t \in I, \min_{r=1, \ldots, j} x_{i_r}(t) \leq y \leq \max_{r=1, \ldots, j} x_{i_r}(t) \}. $$Proporción de Inclusión:
El índice de profundidad basado en bandas para $x$ es:
$$ S_n^{(j)}(x) = \binom{n}{j}^{-1} \sum_{1 \leq i_1 < i_2 < \cdots < i_j \leq n} I\{ G(x) \subseteq B(x_{i_1}, x_{i_2}, \ldots, x_{i_j}) \}, $$
donde $I\{\cdot\}$ es la función indicadora que vale 1 si $G(x)$ (la gráfica de $x$) está contenida en la banda, y 0 en caso contrario. Observe que la cantidad anterior calcula la proporción de bandas que contienen a la curva $x$.Profundidad:
La profundidad de banda para $J \geq 2$ se define como:
$$ S_{n,J}(x) = \sum_{j=2}^{J} S_n^{(j)}(x). $$
Modified Band Depth (MBD)¶
La Modified Band Depth (MBD) extiende la definición de BD al permitir medir qué proporción del intervalo $I$ está contenida dentro de una banda. Esto hace que la métrica sea más robusta a comportamientos locales irregulares.
Conjunto de Inclusión Parcial:
Se define el conjunto donde la función $x$ está dentro de la banda como:
$$ A_j(x; x_{i_1}, \ldots, x_{i_j}) = \{ t \in I : \min_{r=1, \ldots, j} x_{i_r}(t) \leq x(t) \leq \max_{r=1, \ldots, j} x_{i_r}(t) \}. $$Medida Proporcional:
Se mide la proporción de tiempo que $x$ está dentro de la banda:
$$ \lambda_r(A_j(x)) = \frac{\lambda(A_j(x))}{\lambda(I)}, $$
donde $\lambda$ es la medida de Lebesgue.Profundidad Generalizada:
La MBD es:
$$ GS_n^{(j)}(x) = \binom{n}{j}^{-1} \sum_{1 \leq i_1 < \cdots < i_j \leq n} \lambda_r(A_j(x; x_{i_1}, \ldots, x_{i_j})). $$
Finalmente:
$$ GS_{n,J}(x) = \sum_{j=2}^{J} GS_n^{(j)}(x). $$
Implementación del Functional Boxplot¶
Ordenación por Profundidad:
Ordena las curvas de acuerdo con su profundidad BD o MBD, desde la más profunda hasta la menos profunda.Mediana Funcional:
La curva más profunda se define como la mediana funcional:
$$ \hat{m}_n = \arg\max_{x \in \{x_1, \ldots, x_n\}} S_n^{(J)}(x) \quad \text{(o MBD)}. $$Región Central y Outliers:
- Define una banda central basada en el rango intercuartílico de las curvas más profundas.
- Las curvas fuera de esta banda se consideran outliers.
Visualización:
- La mediana funcional y las bandas centrales se visualizan de manera análoga al boxplot estándar, representando las curvas más centrales y las más extremas.
Propiedades y Ventajas¶
Robustez: Ambas métricas son resistentes a valores atípicos, pero la MBD es más flexible frente a irregularidades locales.
Eficiencia Computacional: La BD y MBD son rápidas de calcular, incluso para datos de alta dimensión.
Aplicabilidad: Son útiles para detectar outliers y analizar la centralidad de curvas en conjuntos funcionales.
Implementación: Se encuentra implementada en librerias de R o Python, permitiendo calcular estas métricas y visualizar el functional boxplot.
def agregar_valores(data, values, wavelength):
"""
Agrega un nuevo valor al final de cada fila en la dimensión de wavelength.
Parámetros:
- data: Array numpy con dimensiones (268, 571, ...)
- values: Array numpy con dimensiones (268,)
- wavelength: Índice de la longitud de onda donde se agregará el nuevo valor
Retorna:
- Nuevo array numpy con valores agregados
"""
# Crear un nuevo array con un espacio adicional en la dimensión de wavelength
data_augm = np.zeros((data.shape[0],
data.shape[1] + 1,
data.shape[2]))
# Verificar dimensiones
if values.shape[0] != data.shape[0]:
raise ValueError("Las dimensiones de 'values' y 'data' no coinciden")
# Copiar los datos originales
data_augm[:, :data.shape[1], :] = data
# Agregar los valores al final
for i in range(data.shape[0]):
data_augm[i, -1, wavelength] = values[i]
return data_augm
def analisis_profundidad_funcional(data, fd, wavelength, metodo='mbd'):
"""
Realiza análisis de profundidad funcional y genera visualización de percentiles.
Parámetros:
- data: Array numpy con dimensiones (268, 571, ...)
- fd: Datos funcionales para cálculo de profundidad
- wavelength: Índice de la longitud de onda
- metodo: 'mbd' (ModifiedBandDepth) o 'bd' (BandDepth)
Retorna:
- DataFrame con valores de profundidad ordenados
- Gráfico de percentiles
"""
# Selección del método de profundidad
if metodo == 'mbd':
depth = skfda.exploratory.depth.ModifiedBandDepth()
elif metodo == 'bd':
depth = skfda.exploratory.depth.BandDepth()
else:
raise ValueError("Método de profundidad no válido. Use 'mbd' o 'bd'")
# Cálculo de valores de profundidad
values = depth(fd)
# Agregar valores de profundidad al array original
data_modificado = agregar_valores(data, values, wavelength)
# Crear DataFrame
df = pd.DataFrame(data_modificado[:,:, wavelength])
df.sort_values(by=571, axis=0, ascending=False, inplace=True)
# imprime las 10 primeras filas de df
print(f"Valores de profundidad ordenados - Método {metodo.capitalize()}")
display(df.head(20))
# Cálculo de percentiles
percentile_90 = df.quantile(0.90)[:570]
percentile_95 = df.quantile(0.95)[:570]
# Gráfico de percentiles
plt.figure(figsize=(12, 6))
plt.plot(percentile_90, label='Percentil 90%', color='blue', linestyle='--')
plt.plot(percentile_95, label='Percentil 95%', color='red', linestyle='--')
plt.title(f'Percentiles 90% y 95% - Profundidad {metodo.capitalize()}')
plt.xlabel('Muestras')
plt.ylabel('Profundidad')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
return df
df_bd = analisis_profundidad_funcional(data, fd_analysis, wavelength_idx, metodo='bd')
df_mbd = analisis_profundidad_funcional(data, fd_analysis, wavelength_idx, metodo='mbd')
Valores de profundidad ordenados - Método Bd
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 562 | 563 | 564 | 565 | 566 | 567 | 568 | 569 | 570 | 571 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 26 | 1.906 | 2.055 | 2.294 | 2.658 | 3.165 | 3.816 | 4.624 | 5.627 | 6.872 | 8.397 | ... | 16.591 | 16.396 | 16.286 | 16.302 | 16.418 | 16.553 | 16.668 | 16.818 | 17.069 | 0.046816 |
| 57 | 1.627 | 1.808 | 2.112 | 2.540 | 3.078 | 3.736 | 4.568 | 5.617 | 6.853 | 8.180 | ... | 16.564 | 16.182 | 15.986 | 16.019 | 16.187 | 16.410 | 16.766 | 17.419 | 18.385 | 0.042680 |
| 64 | 1.891 | 1.959 | 2.183 | 2.578 | 3.106 | 3.735 | 4.483 | 5.409 | 6.565 | 7.970 | ... | 16.935 | 16.715 | 16.780 | 17.000 | 17.185 | 17.216 | 17.104 | 16.953 | 16.866 | 0.041897 |
| 123 | 1.660 | 1.886 | 2.191 | 2.551 | 2.971 | 3.490 | 4.156 | 4.997 | 6.023 | 7.219 | ... | 17.802 | 17.634 | 17.488 | 17.346 | 17.277 | 17.372 | 17.643 | 18.008 | 18.363 | 0.039298 |
| 31 | 1.932 | 2.032 | 2.262 | 2.628 | 3.111 | 3.716 | 4.484 | 5.459 | 6.659 | 8.093 | ... | 16.733 | 16.472 | 16.290 | 16.276 | 16.447 | 16.756 | 17.160 | 17.656 | 18.218 | 0.038376 |
| 75 | 1.995 | 2.177 | 2.438 | 2.816 | 3.355 | 4.073 | 4.958 | 5.991 | 7.184 | 8.602 | ... | 16.807 | 16.785 | 16.790 | 16.789 | 16.796 | 16.894 | 17.144 | 17.467 | 17.694 | 0.037677 |
| 36 | 1.940 | 2.146 | 2.480 | 2.939 | 3.506 | 4.160 | 4.918 | 5.853 | 7.059 | 8.554 | ... | 17.275 | 17.381 | 17.550 | 17.657 | 17.623 | 17.511 | 17.501 | 17.740 | 18.218 | 0.036140 |
| 93 | 1.697 | 1.802 | 2.037 | 2.414 | 2.909 | 3.517 | 4.286 | 5.269 | 6.474 | 7.889 | ... | 17.627 | 17.714 | 17.588 | 17.368 | 17.186 | 17.111 | 17.178 | 17.437 | 17.893 | 0.035217 |
| 74 | 1.872 | 1.965 | 2.206 | 2.622 | 3.193 | 3.887 | 4.725 | 5.796 | 7.200 | 8.970 | ... | 16.162 | 16.119 | 16.264 | 16.506 | 16.724 | 16.847 | 16.880 | 16.900 | 16.983 | 0.034714 |
| 59 | 1.746 | 1.968 | 2.299 | 2.721 | 3.209 | 3.767 | 4.440 | 5.291 | 6.367 | 7.691 | ... | 16.393 | 16.341 | 16.243 | 16.101 | 15.934 | 15.793 | 15.775 | 15.966 | 16.352 | 0.034546 |
| 79 | 1.859 | 2.115 | 2.443 | 2.845 | 3.349 | 3.997 | 4.853 | 5.999 | 7.487 | 9.277 | ... | 17.132 | 17.309 | 17.461 | 17.541 | 17.591 | 17.698 | 17.946 | 18.392 | 19.004 | 0.033876 |
| 98 | 1.642 | 1.846 | 2.138 | 2.532 | 3.057 | 3.745 | 4.602 | 5.635 | 6.871 | 8.334 | ... | 16.378 | 16.411 | 16.468 | 16.487 | 16.494 | 16.537 | 16.616 | 16.702 | 16.775 | 0.033876 |
| 32 | 1.648 | 1.797 | 2.060 | 2.470 | 3.030 | 3.702 | 4.453 | 5.299 | 6.303 | 7.528 | ... | 16.816 | 16.854 | 16.888 | 16.875 | 16.862 | 16.961 | 17.251 | 17.731 | 18.323 | 0.033177 |
| 126 | 1.654 | 1.893 | 2.263 | 2.738 | 3.276 | 3.854 | 4.492 | 5.245 | 6.181 | 7.357 | ... | 17.201 | 16.850 | 16.612 | 16.590 | 16.720 | 16.870 | 16.992 | 17.149 | 17.406 | 0.033121 |
| 52 | 1.743 | 1.888 | 2.183 | 2.645 | 3.234 | 3.889 | 4.607 | 5.478 | 6.627 | 8.100 | ... | 16.137 | 16.129 | 16.054 | 15.941 | 15.900 | 15.985 | 16.047 | 15.816 | 15.179 | 0.032730 |
| 35 | 2.002 | 2.110 | 2.362 | 2.819 | 3.474 | 4.247 | 5.078 | 5.987 | 7.047 | 8.330 | ... | 16.389 | 16.181 | 16.098 | 16.166 | 16.403 | 16.771 | 17.128 | 17.328 | 17.353 | 0.032562 |
| 41 | 1.978 | 2.146 | 2.425 | 2.849 | 3.412 | 4.094 | 4.904 | 5.885 | 7.073 | 8.482 | ... | 17.442 | 17.466 | 17.420 | 17.310 | 17.169 | 17.051 | 16.997 | 17.009 | 17.053 | 0.032394 |
| 61 | 1.933 | 2.077 | 2.295 | 2.618 | 3.108 | 3.831 | 4.806 | 5.980 | 7.266 | 8.621 | ... | 17.420 | 17.296 | 17.296 | 17.389 | 17.489 | 17.588 | 17.779 | 18.159 | 18.717 | 0.032282 |
| 34 | 1.958 | 2.099 | 2.325 | 2.683 | 3.209 | 3.897 | 4.719 | 5.690 | 6.880 | 8.334 | ... | 17.422 | 17.298 | 17.194 | 17.179 | 17.250 | 17.391 | 17.624 | 17.962 | 18.362 | 0.032282 |
| 55 | 1.737 | 1.982 | 2.370 | 2.903 | 3.561 | 4.327 | 5.205 | 6.205 | 7.339 | 8.652 | ... | 17.794 | 17.599 | 17.474 | 17.416 | 17.453 | 17.621 | 17.890 | 18.156 | 18.334 | 0.032059 |
20 rows × 572 columns
Valores de profundidad ordenados - Método Mbd
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 562 | 563 | 564 | 565 | 566 | 567 | 568 | 569 | 570 | 571 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 174 | 1.362 | 1.439 | 1.614 | 1.876 | 2.189 | 2.551 | 3.033 | 3.726 | 4.668 | 5.815 | ... | 15.712 | 15.364 | 15.146 | 15.063 | 15.070 | 15.139 | 15.312 | 15.672 | 16.221 | 0.484090 |
| 144 | 1.462 | 1.797 | 1.863 | 2.483 | 2.945 | 3.007 | 3.613 | 4.068 | 4.912 | 6.689 | ... | 15.652 | 16.044 | 16.044 | 16.074 | 17.268 | 17.211 | 15.376 | 15.672 | 17.251 | 0.478247 |
| 78 | 2.024 | 2.117 | 2.277 | 2.538 | 2.941 | 3.516 | 4.269 | 5.212 | 6.386 | 7.833 | ... | 16.231 | 16.252 | 16.280 | 16.320 | 16.392 | 16.519 | 16.751 | 17.144 | 17.677 | 0.477717 |
| 183 | 1.365 | 1.445 | 1.598 | 1.824 | 2.124 | 2.527 | 3.080 | 3.813 | 4.704 | 5.695 | ... | 15.607 | 15.488 | 15.427 | 15.448 | 15.538 | 15.661 | 15.772 | 15.848 | 15.905 | 0.477697 |
| 96 | 2.079 | 2.198 | 2.412 | 2.761 | 3.262 | 3.892 | 4.627 | 5.475 | 6.482 | 7.715 | ... | 16.176 | 16.130 | 16.258 | 16.501 | 16.750 | 16.876 | 16.773 | 16.424 | 15.916 | 0.476309 |
| 43 | 1.734 | 1.904 | 2.129 | 2.426 | 2.825 | 3.367 | 4.093 | 5.035 | 6.202 | 7.575 | ... | 16.413 | 16.335 | 16.276 | 16.220 | 16.188 | 16.240 | 16.438 | 16.804 | 17.301 | 0.475230 |
| 52 | 1.743 | 1.888 | 2.183 | 2.645 | 3.234 | 3.889 | 4.607 | 5.478 | 6.627 | 8.100 | ... | 16.137 | 16.129 | 16.054 | 15.941 | 15.900 | 15.985 | 16.047 | 15.816 | 15.179 | 0.475000 |
| 138 | 1.529 | 1.667 | 1.918 | 2.315 | 2.861 | 3.527 | 4.268 | 5.067 | 5.944 | 6.945 | ... | 16.178 | 16.122 | 16.134 | 16.181 | 16.187 | 16.185 | 16.331 | 16.759 | 17.422 | 0.472879 |
| 27 | 1.824 | 2.006 | 2.313 | 2.709 | 3.163 | 3.705 | 4.405 | 5.328 | 6.501 | 7.911 | ... | 16.060 | 16.091 | 16.072 | 16.111 | 16.263 | 16.425 | 16.424 | 16.198 | 15.850 | 0.470869 |
| 175 | 1.399 | 1.476 | 1.632 | 1.859 | 2.154 | 2.533 | 3.029 | 3.665 | 4.441 | 5.342 | ... | 16.526 | 16.405 | 16.435 | 16.556 | 16.617 | 16.541 | 16.430 | 16.450 | 16.646 | 0.467707 |
| 180 | 1.989 | 1.882 | 1.906 | 2.129 | 2.540 | 3.075 | 3.680 | 4.358 | 5.155 | 6.111 | ... | 16.591 | 16.528 | 16.421 | 16.276 | 16.125 | 16.002 | 15.948 | 15.993 | 16.126 | 0.467657 |
| 65 | 2.115 | 2.065 | 2.121 | 2.358 | 2.797 | 3.392 | 4.093 | 4.912 | 5.915 | 7.143 | ... | 16.268 | 16.199 | 16.073 | 15.990 | 16.019 | 16.114 | 16.116 | 15.866 | 15.351 | 0.467653 |
| 25 | 2.122 | 2.153 | 2.351 | 2.730 | 3.247 | 3.876 | 4.661 | 5.678 | 6.952 | 8.438 | ... | 16.246 | 15.890 | 15.689 | 15.721 | 15.980 | 16.356 | 16.674 | 16.796 | 16.705 | 0.467151 |
| 158 | 1.452 | 1.411 | 1.464 | 1.722 | 2.169 | 3.570 | 2.911 | 3.319 | 5.210 | 6.077 | ... | 16.954 | 16.984 | 17.875 | 15.328 | 16.087 | 16.912 | 15.874 | 16.366 | 17.390 | 0.466007 |
| 22 | 1.836 | 2.036 | 2.395 | 2.882 | 3.454 | 4.125 | 4.959 | 5.988 | 7.185 | 8.527 | ... | 17.130 | 17.183 | 17.113 | 17.008 | 16.953 | 16.958 | 17.018 | 17.161 | 17.398 | 0.464580 |
| 100 | 2.249 | 2.330 | 2.561 | 2.930 | 3.398 | 3.954 | 4.634 | 5.509 | 6.656 | 8.121 | ... | 16.295 | 16.342 | 16.505 | 16.708 | 16.874 | 16.970 | 17.041 | 17.168 | 17.381 | 0.463610 |
| 98 | 1.642 | 1.846 | 2.138 | 2.532 | 3.057 | 3.745 | 4.602 | 5.635 | 6.871 | 8.334 | ... | 16.378 | 16.411 | 16.468 | 16.487 | 16.494 | 16.537 | 16.616 | 16.702 | 16.775 | 0.463585 |
| 181 | 1.269 | 1.508 | 1.783 | 2.075 | 2.399 | 2.791 | 3.273 | 3.853 | 4.547 | 5.406 | ... | 16.897 | 16.397 | 16.099 | 16.144 | 16.434 | 16.714 | 16.705 | 16.237 | 15.372 | 0.462293 |
| 143 | 1.514 | 1.655 | 1.883 | 2.220 | 2.676 | 3.259 | 3.979 | 4.830 | 5.777 | 6.799 | ... | 16.120 | 15.983 | 15.836 | 15.721 | 15.659 | 15.616 | 15.503 | 15.260 | 14.920 | 0.461448 |
| 178 | 1.723 | 1.638 | 1.620 | 1.738 | 2.020 | 2.439 | 2.950 | 3.554 | 4.295 | 5.205 | ... | 16.595 | 16.499 | 16.396 | 16.349 | 16.355 | 16.338 | 16.265 | 16.219 | 16.294 | 0.460649 |
20 rows × 572 columns
Boxplot funcional¶
Construcción de Functional Boxplots¶
Un functional boxplot extiende el concepto del boxplot clásico a datos funcionales y consta de los siguientes elementos:
Región central del 50%: Representa el rango central de las curvas más profundas, delimitada por la banda que contiene el 50% de las curvas más centrales.
$$ C_{0.5} = \{ (t, y(t)) : \min_{r=1,\ldots,\lfloor n/2 \rfloor} y_{[r]}(t) \leq y(t) \leq \max_{r=1,\ldots,\lfloor n/2 \rfloor} y_{[r]}(t) \}. $$Curva mediana funcional: La curva más central basada en la profundidad (por ejemplo, BD o MBD).
$$ \hat{m}_n = \arg\max_{y \in \{y_1, \ldots, y_n\}} \text{BD}_n(y). $$Región máxima no atípica: El rango completo de curvas no consideradas outliers.
Valores atípicos: Curvas fuera de los límites establecidos por la región central ampliada en 1.5 veces su rango.
Procedimiento de Construcción¶
Ordenación: Las curvas se ordenan según su profundidad (BD o MBD) de forma descendente.
Cálculo de la región central: Se identifica la banda que contiene las curvas más profundas (usualmente el 50%).
Detección de outliers: Se identifican curvas fuera de los límites ampliados.
Visualización:
- La región central del 50% se representa como el análogo funcional del rango intercuartílico (IQR).
- La mediana funcional se muestra como la curva más representativa.
- Los valores atípicos se destacan.
Conclusiones¶
- Robustez: El functional boxplot es robusto frente a valores atípicos y proporciona una representación visual confiable de la dispersión y centralidad de las curvas.
- Descriptivo y Comparativo: Permite comparar poblaciones de curvas y detectar diferencias significativas entre conjuntos de datos funcionales.
- Aplicaciones Versátiles: Útil en múltiples campos como datos temporales, espaciales y espaciotemporales.
def functional_boxplot(data, wavelength_idx, F=1.5):
"""
Create functional boxplot with outliers
"""
X = data[:, :, wavelength_idx]
fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio, dataset_name=f'Boxplot with F: {F}')
fdBoxplot = Boxplot(fd, factor=F)
fdBoxplot.plot()
return fdBoxplot
fboxplot = functional_boxplot(data, wavelength_idx)
outliers = fboxplot.outliers
# Compara los outliers con data y muestrame los valores de data que son outliers
data_outliers = pd.DataFrame(data[outliers, :, wavelength_idx])
display(data_outliers)
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 561 | 562 | 563 | 564 | 565 | 566 | 567 | 568 | 569 | 570 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1.891 | 2.045 | 2.286 | 2.603 | 2.995 | 3.497 | 4.161 | 5.039 | 6.181 | 7.640 | ... | 25.238 | 25.418 | 25.498 | 25.481 | 25.450 | 25.457 | 25.492 | 25.541 | 25.645 | 25.852 |
| 1 | 1.966 | 2.171 | 2.497 | 2.931 | 3.456 | 4.090 | 4.863 | 5.796 | 6.910 | 8.263 | ... | 22.321 | 22.459 | 22.633 | 22.694 | 22.543 | 22.194 | 21.763 | 21.392 | 21.164 | 21.059 |
| 2 | 1.206 | 1.321 | 1.466 | 1.650 | 1.884 | 2.183 | 2.580 | 3.117 | 3.816 | 4.655 | ... | 24.670 | 24.268 | 24.006 | 23.916 | 23.959 | 24.051 | 24.130 | 24.230 | 24.465 | 24.887 |
3 rows × 571 columns
Boxplot funcional adjustado¶
Ajuste del Factor en Functional Boxplots¶
El factor constante 1.5, utilizado tradicionalmente en boxplots clásicos para la detección de outliers, también se aplica en functional boxplots como multiplicador del rango de la región central del 50%. Sin embargo, este valor fijo puede no ser adecuado en todos los contextos debido a la naturaleza de los datos funcionales.
Razones para Ajustar el Factor¶
Dependencia Temporal y Espaciotemporal:
- Las curvas funcionales suelen presentar dependencia temporal (en datos univariados) y espaciotemporal (en datos multivariados).
- Esta dependencia puede alterar la dispersión observada en las curvas y, por tanto, afecta la detección de outliers si se usa un factor fijo como 1.5.
Distribución de los Outliers:
- En un boxplot clásico, el factor 1.5 está basado en la distribución normal estándar, donde la probabilidad de detectar un valor como outlier es aproximadamente 0.7%.
- Este factor no toma en cuenta la correlación inherente en los datos funcionales, lo que puede llevar a sobreestimación o subestimación de outliers.
Impacto en la Visualización y Detección:
- Un factor fijo puede no ser adecuado si las curvas tienen una estructura de correlación compleja.
- Ajustar el factor permite que la probabilidad de detectar valores atípicos se adapte mejor a las características del conjunto de datos.
Propuesta de Ajuste¶
El factor puede ajustarse dinámicamente para garantizar que la probabilidad de no detectar outliers sea consistente con la estructura de los datos. Por ejemplo:
- En ausencia de outliers reales, ajustar el factor para que la probabilidad de no detectar outliers sea del 99.3%.
- Este enfoque mejora la sensibilidad del functional boxplot para identificar valores atípicos mientras mantiene la robustez frente a la dependencia.
Simulación¶
Simularemos curvas, a apartir de un proceso Gaussiano, definiendo la estructura de variación a partir de una matriz de correlaciones basadas en un estimador robusto de la varianza y covarianza dado por la siguiente espresión:
Estimador robusto de una matriz de dispersión (Ma & Genton, 2001)¶
El estimador de covarianza robusta está dado por:
$$ \text{Cov}(X, Y) = \frac{\sigma_X \sigma_Y}{4} \left[ \text{Var}\left( \frac{X}{\sigma_X} + \frac{Y}{\sigma_Y} \right) - \text{Var}\left( \frac{X}{\sigma_X} - \frac{Y}{\sigma_Y} \right) \right], $$
Este estimador reduce la influencia de outliers en la matriz de dispersión, crucial en escenarios con dependencia estructural en los datos.
Estimador altamente robusto de escala (Rousseeuw y Croux, 1992, 1993)¶
Dado un conjunto de datos $\{z_1, \ldots, z_n\}$, el estimador robusto de escala es:
$$ \hat{\sigma}_Z = Q_n(Z) = d \, \Big[ \lvert z_i - z_j \rvert; \, i < j, \, i, j = 1, 2, \ldots, n \Big]_{(k)}, $$
donde:
- $d = 2.2191$ es un factor de consistencia,
- $k = \lfloor ((n/2) + 2)/4 \rfloor + 1$,
- $(k)$ indica aproximadamente el primer cuartil para muestras grandes.
Este método permite medir la escala de los datos mientras mantiene robustez frente a valores atípicos.
Una vez que se han calculado la matriz $Cov(X,Y)$ con la data inicial, se procede a hacer una serie de simulaciones para encontrar el mejor F. Para esto simulamos 1000 matrices de muestras de los datos, y calculamos
Conclusión¶
El ajuste del factor en functional boxplots es crucial para asegurar una detección de outliers confiable en escenarios con alta dependencia temporal o espaciotemporal. Aunque el valor 1.5 funciona bien en contextos simples, en datos funcionales complejos, un factor adaptable proporciona mejores resultados sin comprometer la visualización.
import numpy as np
import skfda
from skfda.exploratory.outliers import BoxplotOutlierDetector
import numpy as np
import scipy.special as sc
def calculate_stats(data):
"""
Calcula un cuantil específico de las diferencias absolutas para cada columna
de una matriz 2D, considerando únicamente los pares donde i < j.
Args:
data (numpy.ndarray): Una matriz 2D donde cada columna representa una variable.
Returns:
numpy.ndarray: Un arreglo con los cuantiles calculados para cada columna.
Raises:
ValueError: Si la entrada no es una matriz 2D numérica.
"""
if not isinstance(data, np.ndarray) or data.ndim != 2:
raise ValueError("Input data must be a 2D NumPy array.")
if not np.issubdtype(data.dtype, np.number):
raise ValueError("Input data must be numeric.")
n_cols = data.shape[1]
d = 2.2191
quantiles = []
for col_idx in range(n_cols):
col_data = data[:, col_idx]
n = len(col_data)
diffs = []
# Calcular diferencias solo para i < j
for i in range(n):
for j in range(i + 1, n): # Solo i < j
diffs.append(np.abs(col_data[i] - col_data[j]))
diffs = np.array(diffs)
# Calcular el índice del cuantil específico
k = int(np.floor((sc.comb(n, 2) + 2) / 4)) + 1
# multiplica k por d
k_int = int(k * d)
quantiles.append(np.partition(diffs, k - 1)[k - 1]) # Obtener el k-ésimo valor más pequeño
return np.array(quantiles)
def calculate_covariance_matrix(array_reshaped, stats):
# Número de columnas
num_columns = array_reshaped.shape[1]
# Inicializar la matriz de covarianza
covariance_matrix = np.zeros((num_columns, num_columns))
for i in range(num_columns):
for j in range(num_columns):
# Extraer las columnas i y j
X = array_reshaped[:, i, 4]
Y = array_reshaped[:, j, 4]
# Obtener sigma_X y sigma_Y de stats
sigma_X = stats[i]
sigma_Y = stats[j]
# Calcular términos de la fórmula
term_1 = (X / sigma_X) + (Y / sigma_Y)
term_2 = (X / sigma_X) - (Y / sigma_Y)
# Calcular las varianzas
var_term_1 = np.var(term_1)
var_term_2 = np.var(term_2)
# Aplicar la fórmula
covariance_matrix[i, j] = (sigma_X * sigma_Y / 4) * (var_term_1 - var_term_2)
return covariance_matrix
def calculate_outlier_proportions(covariance_matrix, num_matrices=1000, num_samples=100, num_groups=10):
# Generar las 1000 matrices de muestras
concatenated_samples = [
np.random.multivariate_normal(np.zeros(covariance_matrix.shape[0]), covariance_matrix, size=num_samples)
for _ in range(num_matrices)
]
# Dividir las matrices en 10 grupos
matrices_per_group = num_matrices // num_groups
groups = [
concatenated_samples[i * matrices_per_group:(i + 1) * matrices_per_group]
for i in range(num_groups)
]
# Configurar los valores de factor para cada grupo
factors = [0.9 + 0.1 * i for i in range(num_groups)]
# Inicializar resultados
group_proportions = []
# Procesar cada grupo de matrices
for group_idx, (group, factor) in enumerate(zip(groups, factors)):
proportions = []
for sample in group:
# Crear el objeto FDataGrid para la matriz actual
num_columns = covariance_matrix.shape[0]
grid_points = np.arange(275, 275 + 0.5 * num_columns, 0.5)
fd = skfda.FDataGrid(data_matrix=sample, grid_points=grid_points)
# Aplicar el detector de outliers
out_detector = BoxplotOutlierDetector(factor=factor)
outlier_predictions = out_detector.fit_predict(fd)
# Calcular proporción de outliers
proportion = np.sum(outlier_predictions == -1) / len(outlier_predictions)
proportions.append(proportion)
# Calcular promedio de proporciones para el grupo actual
group_avg_proportion = 1-np.mean(proportions)
group_proportions.append(group_avg_proportion)
print(f"Grupo {group_idx + 1} con factor {factor:.1f}: Promedio de proporciones de outliers = {group_avg_proportion:.4f}")
return group_proportions
stats = calculate_stats(data[:,:,wavelength_idx])
covariance_matrix = calculate_covariance_matrix(data, stats)
calculate_outlier_proportions(covariance_matrix, num_matrices=1000, num_samples=100, num_groups=10)
Grupo 1 con factor 0.9: Promedio de proporciones de outliers = 0.9603 Grupo 2 con factor 1.0: Promedio de proporciones de outliers = 0.9659 Grupo 3 con factor 1.1: Promedio de proporciones de outliers = 0.9782 Grupo 4 con factor 1.2: Promedio de proporciones de outliers = 0.9878 Grupo 5 con factor 1.3: Promedio de proporciones de outliers = 0.9906 Grupo 6 con factor 1.4: Promedio de proporciones de outliers = 0.9941 Grupo 7 con factor 1.5: Promedio de proporciones de outliers = 0.9976 Grupo 8 con factor 1.6: Promedio de proporciones de outliers = 0.9983 Grupo 9 con factor 1.7: Promedio de proporciones de outliers = 0.9989 Grupo 10 con factor 1.8: Promedio de proporciones de outliers = 0.9989
[np.float64(0.9603), np.float64(0.9659), np.float64(0.9782), np.float64(0.9878), np.float64(0.9906), np.float64(0.9941), np.float64(0.9976), np.float64(0.9983), np.float64(0.9989), np.float64(0.9989)]
fboxplot = functional_boxplot(data, wavelength_idx, F=1.4)
outliers = fboxplot.outliers
# Compara los outliers con data y muestrame los valores de data que son outliers
data_outliers = pd.DataFrame(data[outliers, :, wavelength_idx])
Análisis de la Profundidad de Variación Total y Similitud de Forma¶
** 1. Introducción**¶
Sea $X$ un proceso estocástico en $\tau$, donde $\tau$ es un intervalo en $\mathbb{R}$, con distribución $F_X$. Denotamos:
- $f$: una función.
- $f(t)$: el valor funcional en un tiempo $t$ dado.
2. Profundidad de Variación Total Puntual¶
Para una función $f(t)$, definimos:
$$ R_f(t) = \mathbb{I}_{\{ X(t) \leq f(t) \}}, $$
donde $\mathbb{I}$ es la función indicadora. A partir de esta relación:
$$ p_f(t) = \mathbb{E}[R_f(t)] = \mathbb{P}\{ X(t) \leq f(t) \}, $$
lo que indica la posición de $f(t)$ con respecto a $X(t)$.
La profundidad de variación total puntual se define como:
$$ D_f(t) = \text{Var}[R_f(t)] = p_f(t)(1 - p_f(t)). $$
3. Descomposición de la Profundidad de Variación Total¶
Sea $s, t$ dos puntos en el tiempo donde $s = t - \Delta$. La profundidad de variación total puntual tiene la siguiente descomposición:
$$ D_f(t) = \text{Var}[R_f(t)] = \text{Var}[\mathbb{E}\{R_f(t) \mid R_f(s)\}] + \mathbb{E}[\text{Var}\{R_f(t) \mid R_f(s)\}]. $$
Esta descomposición implica que la varianza total de $R_f(t)$ se puede dividir en:
- $\text{Var}[\mathbb{E}\{R_f(t) \mid R_f(s)\}]$: la componente de forma, que representa la porción de variabilidad de $R_f(t)$ explicada por $R_f(s)$.
- $\mathbb{E}[\text{Var}\{R_f(t) \mid R_f(s)\}]$: la componente de magnitud, que es independiente de $R_f(s)$.
4. Profundidad de Variación Total Funcional (TVD)¶
La profundidad de variación total funcional (TVD, por sus siglas en inglés) para una función $f$ es:
$$ TVD(f) = \int w(t) D_f(t) \, dt, $$
donde $w(t)$ es una función de peso definida sobre el dominio de $f(t)$. Si consideramos que $w(t)$ es constante, tal que:
$$ w(t) \equiv \frac{1}{|\tau|}, $$
entonces el TVD se comporta de manera similar a la profundidad de banda modificada.
5. Similitud de Forma (Shape Similarity, SS)¶
Para cualquier función $f$ en $\mathbb{R}$, la similitud de forma dada una ventana temporal $\Delta$ se define como:
$$ SS(f; \Delta) = \int v(t; \Delta) S_f(t; \Delta) \, dt, $$
donde:
- $v(t; \Delta)$: función de peso, dada por:
$$ v(t; \Delta) = \frac{|f(t) - f(t - \Delta)|}{\int |f(t) - f(t - \Delta)|}. $$
- $S_f(t; \Delta)$: componente de forma, definida como:
$$ S_f(t; \Delta) = \begin{cases} \frac{\text{Var}[R_f(t) \mid R_f(t - \Delta)]}{D_f(t)}, & \text{si } D_f(t) \neq 0, \\ 1, & \text{si } D_f(t) = 0. \end{cases} $$
Interpretación de la Similitud de Forma¶
Valores pequeños de $SS(f; \Delta)$ están asociados con una mayor desviación en la forma. Sin embargo, para pares atípicos $(f(t - \Delta), f(t))$, si el denominador $D_f(t)$ es muy pequeño, $S_f(t; \Delta)$ podría no reflejar adecuadamente la desviación.
Para capturar mejor esta desviación, se centra el par $(f(t - \Delta), f(t))$ mediante un desplazamiento $\delta_t$, donde:
$$ \delta_t = f(t) - \text{mediana}\{X(t)\}. $$
El par ajustado es:
$$ (f(t - \Delta), f(t)) - \delta_t. $$
6. Similitud de Forma Modificada (Modified Shape Similarity, MSS)¶
La similitud de forma modificada (MSS) para una función $f$, dada una ventana temporal fija $\Delta$, se define como:
$$ MSS(f; \Delta) = \int v(t; \Delta) S_f(t; \Delta) \, dt, $$
donde $S_f(t; \Delta)$ es:
$$ S_f(t; \Delta) = \frac{\text{Var}[R_f(t) \mid R_f(t - \Delta)]}{D_f(t)}. $$
Además, se define $f(s; \Delta)$ como:
$$ f(s; \Delta) = \begin{cases} f(s) - f(s + \Delta) + \text{mediana}\{X(s + \Delta)\}, & \text{si } s \in \mathbb{R}. \end{cases} $$
Outlier Detection Procedure¶
The process for detecting outliers among a set of $n$ sample curves involves the following steps:
Compute Metrics
- Estimate the total variation depth and modified shape similarity for each curve as before
Identify Shape Outliers
- Create a classical boxplot for the $n$ values of the modified shape similarity.
- Detect outliers using the $F \times IQR$ empirical rule, where $F$ is a user-adjustable factor.
- Curves with modified shape similarity values below the lower fence are classified as shape outliers.
Identify Magnitude Outliers
- Remove shape outliers and construct a functional boxplot using the total variation depth.
- Detect magnitude outliers by identifying curves outside the 50% central region (based on the original observations) expanded by a factor of 1.5. This factor can be adjusted using bootstrap methods .
import rpy2
%load_ext rpy2.ipython
X = data[:, :, wavelength_idx]
%R -i X
Además: Aviso: In (function (package, help, pos = 2, lib.loc = NULL, character.only = FALSE, : library ‘/usr/lib/R/site-library’ contains no packages
%%R
library("fdaoutlier")
# Just show me the first 5 rows in R
class_outliners <- tvdmss(X)
%R -o class_outliners
def plot_kind_class_outliners(outliers, X):
"""
Plot outliers and non-outliers with different colors and labels
"""
# Convierte los índices de outliers a un conjunto para una búsqueda eficiente
magnitude_outliers_set = set(outliers['magnitude_outliers'])
shape_outliers_set = set(outliers['shape_outliers'])
# Crear máscaras para diferentes tipos de outliers
magnitude_outliers_mask = np.array([i in magnitude_outliers_set for i in range(len(X))])
shape_outliers_mask = np.array([i in shape_outliers_set for i in range(len(X))])
# Encuentra los magnitude outliers en X
magnitude_outliers = X[magnitude_outliers_mask, :]
# Encuentra los shape outliers en X
shape_outliers = X[shape_outliers_mask, :]
# Encuentra las curvas que no son outliers
non_outliers_mask = ~(magnitude_outliers_mask | shape_outliers_mask)
non_outliers = X[non_outliers_mask, :]
# Grafica los tipos de curvas
plt.figure(figsize=(12, 6))
for i in range(shape_outliers.shape[0]):
plt.plot(shape_outliers[i], color='blue', alpha=0.7, label=('' if i==0 else '_') + 'Outliers de Forma')
for i in range(magnitude_outliers.shape[0]):
plt.plot(magnitude_outliers[i], color='red', alpha=0.7, label=('' if i==0 else '_') + 'Outliers de Magnitud')
plt.title('Outliers de Magnitud y Forma')
plt.xlabel('Longitud de onda (nm)')
plt.ylabel('Intensidad')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# Llamada a la función
plot_kind_class_outliners(class_outliners, X)
# imprime cuantos outliers de magnitud y forma hay
display(class_outliners['magnitude_outliers'])
print(f"Outliers de magnitud: {len(class_outliners['magnitude_outliers'])}")
print(f"Outliers de forma: {len(class_outliners['shape_outliers'])}")
array([ 10, 71, 131], dtype=int32)
Outliers de magnitud: 3 Outliers de forma: 22
Reproducir figura 2 del paper
Mediana funcional multivariada¶
Definición de la Profundidad Funcional Multivariada (MFD)¶
Proceso Considerado:
- Se analiza un proceso estocástico $K$-variado $\{Y(t), t \in U\}$ definido en $\mathbb{R}^K$, con función de distribución acumulada $F_Y$ que genera caminos continuos en $C(U)^K$ (espacio de funciones continuas).
Medida de Profundidad:
- La MFD para una curva arbitraria $X \in C(U)^K$ se define como: $$ \text{MFD}(X; F_Y) = \int_U D(X(t); F_{Y(t)}) \cdot w(t) \, dt, $$ donde:
- $D$: Función de profundidad estadística en $\mathbb{R}^K$.
- $w(t)$: Función de peso definida en $U$, que se normaliza a $1$, es decir, $\int_U w(t) dt = 1$.
Opciones para la Función de Peso:
Regiones de Profundidad¶
- Definición: La región de profundidad $ D_\alpha(F_X) $ en el nivel $ \alpha \geq 0 $ se define como: $$ D_\alpha(F_X) = \{ x \in \mathbb{R}^K : D(x; F_X) \geq \alpha \}, $$ donde $ D(x; F_X) $ es una función de profundidad que mide cuán central es $ x $ respecto a la distribución $ F_X $. Valores mayores de $ D(x; F_X) $ indican mayor centralidad.
Funciones de Peso¶
Las funciones de peso $ w(t) $ ajustan la contribución de diferentes puntos en el tiempo $ t $ en el análisis de datos funcionales. Se presentan dos ejemplos:
Peso Constante con Función Indicadora:
- Un peso constante dentro de un rango de interés permite excluir ciertos periodos, como:
- Fases de arranque en procesos industriales.
- Mediciones imprecisas en regiones específicas de datos espectrales.
- Un peso constante dentro de un rango de interés permite excluir ciertos periodos, como:
Peso Ajustado por Variabilidad:
- Una función de peso que se adapta a la variabilidad local de los datos: $$ w(t) = w_\alpha(t; F_Y(t)) = \frac{\text{vol}(D_\alpha(F_Y(t)))}{\int_U\text{vol}(D_\alpha(F_{Y(u)})) \, du} $$ donde $ \text{vol}(D_\alpha(F_Y(t))) $ mide el "volumen" de la región de profundidad $ D_\alpha(F_Y(t)) $ en el tiempo $ t $. Este enfoque considera la variabilidad en amplitud (variabilidad vertical) de los datos funcionales.
Definición de Funciones y Correspondencias¶
Función $H$:
- Se define como $H : U \times \mathbb{R}^K \to [0, 1]$: $$ H(t, x) = D(x, F_Y(t)), $$ donde $D(x, F_Y(t))$ es la profundidad del vector $x$ respecto a la distribución $F_Y(t)$.
Correspondencia $\gamma$:
- Para cada $t \in U$, $\gamma(t)$ es el conjunto de $x \in \mathbb{R}^K$ en los que se evalúa la profundidad. Este conjunto es compacto y no vacío.
Funciones auxiliares:
- $\Gamma: U \to \mathbb{R}$ se define como: $$ \Gamma(t) = (-\infty, \max_{x \in \gamma(t)} D(x; F_Y(t))]. $$
- $G(t)$ denota el gráfico de $H(t, \cdot)$ restringido a $\gamma(t)$, y su clausura en $\mathbb{R}^K \times \mathbb{R} \cup \{-\infty\}$ se denota como $G(t)$.
Teorema 2: Propiedades del Máximo de Profundidad (Existencia de la mediana)¶
Suposiciones:
- $H$ es semicontinua superior.
- $\gamma$ es una correspondencia compacta y semicontinua superior.
- $\Gamma$ es semicontinua inferior.
Conclusiones:
- La función $t \mapsto \max_{x \in \gamma(t)} D(x; F_Y(t))$ es continua.
- El conjunto de puntos donde la profundidad máxima se alcanza: $$ \Pi(t) = \text{argmax}_{x \in \gamma(t)} D(x; F_Y(t)) $$ es compacto, no vacío y semicontinuo superior. Si $\Pi(t)$ es unívoco, es continuo.
- Para una curva $\vartheta(t) \in \Pi(t)$ y $w(t)$ una función de peso, se tiene: $$ \text{MFD}(X; F_Y) = \int_U D(\vartheta(t); F_Y(t)) \cdot w(t) \, dt. $$
- Si $\vartheta \in C(U)^K$ y maximiza la MFD, entonces $\vartheta(t) \in \Pi(t)$ en cada $t$.
La profundidad resultante se denomina profundidad de semiespacio funcional multivariante (MFHD).
Definición de Muestras Finitas¶
Para observaciones de curvas multivariantes $\{Y^1(t_j), \dots, Y^N(t_j)\}$ en tiempos discretos $t_1, \dots, t_T$:
Profundidad Muestral Multivariante Funcional (MFD): $$ \text{MFD}_N(X) = \sum_{j=1}^T D(X(t_j); F_Y(t_j), N) \cdot W_j, $$ donde $W_j$ es un peso que depende de $t_j$ y $w(t)$.
Consistencia: Se demuestra que $\text{MFD}_N$ converge a $\text{MFD}$ cuando tanto $N$ como $T$ tienden a infinito.
La media se define como la curva más profunda que para MFD por teorema 2 existe
lambdas_3 = data[:, :, 4:]
display(lambdas_3.shape)
# Cambia la dimension de lambdas_3 intercambia por (571, 268, 3)
lambdas_3_new = np.moveaxis(lambdas_3, 0, 1)
display(lambdas_3_new.shape)
%R -i lambdas_3_new
(268, 571, 3)
(571, 268, 3)
%%R
library("mrfDepth")
median_mf <- mrfDepth::mfdmedian(lambdas_3_new)
median_mf
$MFDmedian
[,1] [,2] [,3]
[1,] 2.105005 0.8049057 1.587502
[2,] 2.187033 0.8262345 1.696407
[3,] 2.337496 0.8500005 1.888822
[4,] 2.543371 0.8977931 2.194584
[5,] 2.808741 0.9457164 2.610831
[6,] 3.109565 1.0069587 3.140136
[7,] 3.441962 1.0796899 3.748493
[8,] 3.791654 1.1650708 4.530298
[9,] 4.187731 1.2702830 5.472409
[10,] 4.602064 1.4289436 6.637180
[11,] 5.038211 1.5979974 7.921381
[12,] 5.506924 1.8293113 9.483342
[13,] 5.986744 2.0747602 11.140498
[14,] 6.473741 2.3481041 13.186854
[15,] 6.968772 2.6563346 15.320844
[16,] 7.523273 2.9698299 17.749692
[17,] 8.128189 3.3621962 20.502031
[18,] 8.711750 3.8187841 23.694310
[19,] 9.354988 4.2993876 26.936582
[20,] 9.998370 4.7912929 30.372750
[21,] 10.604316 5.4053845 34.377443
[22,] 11.369136 6.0161060 38.780952
[23,] 12.138011 6.7123730 43.627747
[24,] 12.886487 7.4217204 48.598987
[25,] 13.738171 8.2579607 53.902744
[26,] 14.575708 9.1407035 60.028082
[27,] 15.450214 10.0540269 66.105291
[28,] 16.374814 11.0654319 72.785523
[29,] 17.244564 11.9764134 78.900961
[30,] 18.312495 13.2166005 86.890952
[31,] 19.022477 14.1664380 92.618006
[32,] 20.116644 15.3374546 100.484675
[33,] 20.946326 16.3419502 107.213824
[34,] 22.117082 17.6666607 115.869443
[35,] 23.123445 18.8453012 121.897584
[36,] 24.257791 20.1597170 131.805563
[37,] 25.279525 21.3354233 137.616299
[38,] 26.525818 22.7435449 147.237739
[39,] 27.481334 24.0629645 153.810175
[40,] 28.934868 25.5482467 163.883248
[41,] 30.132987 27.1378440 170.173566
[42,] 31.489360 28.5480889 179.655211
[43,] 32.566663 29.6490323 187.907950
[44,] 33.689635 31.0696740 195.343560
[45,] 35.002795 32.5226736 204.341586
[46,] 36.070027 33.8932282 211.257541
[47,] 37.405291 35.3533386 219.636274
[48,] 38.315319 36.6344593 225.285407
[49,] 39.733432 38.0722759 231.668833
[50,] 40.716470 39.4451926 237.818857
[51,] 41.880403 40.8178573 243.546478
[52,] 42.755958 42.0456595 249.612161
[53,] 43.729239 42.9736918 254.974255
[54,] 44.832729 44.1988506 260.858615
[55,] 45.769214 45.2852136 265.970954
[56,] 46.692657 46.4081007 271.443775
[57,] 47.645147 47.5497123 274.399993
[58,] 48.479773 48.5803963 279.004629
[59,] 49.302679 49.5993938 282.368861
[60,] 50.023688 50.5665519 284.815622
[61,] 50.642496 51.4793485 287.497801
[62,] 51.385963 52.3781814 289.308648
[63,] 52.001433 53.0314963 291.130580
[64,] 52.465882 53.7300048 292.176310
[65,] 53.070219 54.4448569 293.609875
[66,] 53.414099 55.1985523 295.837189
[67,] 54.230735 56.1210323 296.862386
[68,] 54.715249 56.9990554 298.768200
[69,] 55.508267 57.8469431 299.003691
[70,] 56.411754 58.4437657 300.650862
[71,] 57.126491 59.2600650 300.777101
[72,] 57.734251 60.0987221 300.711123
[73,] 58.382562 60.8723601 301.878735
[74,] 59.104979 61.6448273 301.763763
[75,] 59.765328 62.3160232 301.530338
[76,] 60.864386 63.0506014 301.295772
[77,] 61.668376 63.7861976 301.288113
[78,] 62.464861 64.5770600 300.372193
[79,] 63.372372 65.2212285 299.885478
[80,] 64.371142 65.8720298 299.164187
[81,] 65.209276 66.8523518 297.668518
[82,] 66.444708 67.5971799 296.686262
[83,] 67.169503 67.9478712 294.661078
[84,] 68.190111 68.8649447 294.570672
[85,] 68.939222 69.4609354 291.548062
[86,] 70.113131 70.1205653 290.914961
[87,] 71.308014 70.9753590 291.068861
[88,] 72.046522 71.6525559 289.301283
[89,] 73.522340 72.4178235 289.256922
[90,] 74.553856 73.4023682 288.420932
[91,] 75.677896 74.3414597 288.487399
[92,] 77.170862 75.0079096 287.501528
[93,] 78.772393 76.0035190 285.077186
[94,] 79.837337 76.5418750 285.063474
[95,] 81.120451 77.3082635 283.418505
[96,] 82.622196 78.1044900 283.894074
[97,] 83.995226 78.9419488 281.893248
[98,] 85.160245 79.6803612 280.786553
[99,] 86.268862 80.5397237 279.399074
[100,] 87.302130 81.4258515 279.237465
[101,] 88.718670 82.0262374 277.435263
[102,] 90.015457 82.9256042 276.858727
[103,] 91.442713 83.8977394 275.422417
[104,] 92.697041 84.5633156 275.782874
[105,] 94.202897 85.4767722 273.763984
[106,] 95.606896 86.3002161 274.466289
[107,] 97.050653 87.0737461 273.813269
[108,] 98.432641 88.0843678 272.783111
[109,] 99.549037 88.7392417 271.490510
[110,] 100.975330 89.4130086 270.759227
[111,] 102.340477 90.1847089 269.369980
[112,] 103.658406 90.6050729 267.694011
[113,] 104.772167 91.1457487 266.702638
[114,] 106.193843 91.8482942 265.076138
[115,] 107.112255 92.3406410 263.160936
[116,] 108.356982 92.8690691 262.682393
[117,] 109.393185 93.5200331 260.829866
[118,] 110.582410 94.0708591 259.663137
[119,] 111.728098 94.6709009 259.688548
[120,] 112.860839 95.2643751 259.362386
[121,] 113.708150 95.4860430 258.483040
[122,] 115.293849 96.2248730 256.820611
[123,] 115.781633 96.5562556 256.593981
[124,] 116.458086 96.9907500 254.959935
[125,] 117.323280 97.4529138 253.843093
[126,] 118.531773 97.8266076 252.890896
[127,] 119.570154 98.3498613 252.299779
[128,] 120.367064 98.6225334 251.138060
[129,] 121.241266 99.1094153 250.870792
[130,] 121.835021 99.3150120 249.447604
[131,] 122.697827 99.5929174 249.070554
[132,] 122.907830 99.5181049 247.592968
[133,] 123.691463 99.9772003 247.127022
[134,] 124.531578 100.3340044 246.771631
[135,] 124.865036 100.4246486 245.558552
[136,] 125.334134 100.6581817 244.648907
[137,] 125.717981 100.6112761 243.923803
[138,] 126.270852 100.9430542 242.730814
[139,] 126.633448 101.1598014 242.061484
[140,] 126.976057 101.3887108 241.154035
[141,] 127.450503 101.5197168 240.537018
[142,] 127.950741 101.4732367 239.801896
[143,] 128.372640 101.8668673 239.261676
[144,] 128.682900 101.8953378 238.077881
[145,] 128.826146 101.8710395 237.458382
[146,] 128.740573 101.9278875 236.270508
[147,] 129.096224 102.1454971 235.586336
[148,] 129.305024 102.1706923 234.748144
[149,] 129.615068 102.3706690 234.024475
[150,] 129.829666 102.3768041 233.263806
[151,] 129.965066 102.4465602 232.525054
[152,] 130.028368 102.5625015 231.895422
[153,] 130.184601 102.6410795 231.059602
[154,] 130.318059 102.7343891 230.503040
[155,] 130.956052 103.0710442 230.224164
[156,] 130.954831 103.1024725 229.469580
[157,] 130.855169 103.1192938 228.568532
[158,] 130.622216 102.9142646 227.575026
[159,] 130.862398 103.2085299 227.260008
[160,] 130.835379 103.1736036 226.545650
[161,] 130.664783 102.9567262 225.752559
[162,] 130.312057 102.8380172 224.990103
[163,] 130.300527 102.8519813 224.516414
[164,] 130.521593 103.3841075 224.303310
[165,] 131.333597 104.0278120 224.287858
[166,] 130.585639 103.4720112 222.829166
[167,] 130.418296 103.6329203 222.137526
[168,] 129.413526 103.2753928 222.514421
[169,] 129.045711 103.1047861 221.136262
[170,] 128.726833 103.0796486 219.958045
[171,] 128.803166 103.1895023 218.911147
[172,] 128.233996 103.1147451 218.073913
[173,] 127.693202 103.0655587 217.130998
[174,] 127.417675 103.0134163 215.526439
[175,] 126.835367 102.6298819 215.494891
[176,] 125.881504 102.6107986 214.936164
[177,] 125.322834 102.5139776 214.045149
[178,] 125.002858 102.3793994 212.235979
[179,] 124.742566 102.2809006 210.825262
[180,] 124.409985 102.2239233 209.700114
[181,] 123.417046 102.1853923 209.162352
[182,] 122.674101 101.9789017 208.718170
[183,] 121.820939 101.8772595 207.963871
[184,] 121.724441 101.6641297 206.195380
[185,] 121.020380 101.5491043 205.097799
[186,] 120.910624 101.6130157 203.361406
[187,] 120.276186 101.2625467 202.059486
[188,] 119.370869 101.1271802 201.569717
[189,] 119.063494 100.9183789 200.547428
[190,] 118.167734 100.7582471 198.732096
[191,] 116.792695 100.5839388 198.949333
[192,] 116.497655 100.5694041 198.332284
[193,] 115.527783 100.6114467 196.932517
[194,] 115.024928 100.4023104 195.179109
[195,] 114.417714 100.2613160 194.731465
[196,] 113.909470 100.2736321 193.710491
[197,] 113.434083 99.9753197 192.559474
[198,] 112.821986 99.8383121 190.763013
[199,] 111.537973 99.7128760 189.657999
[200,] 110.289522 99.5556736 189.331370
[201,] 110.077287 99.6041187 187.443364
[202,] 108.887043 99.3398761 187.363985
[203,] 109.207405 99.3477490 184.903670
[204,] 108.664330 98.9232762 183.686359
[205,] 107.729279 98.8654745 183.048823
[206,] 106.711420 98.6172655 181.843206
[207,] 105.970114 98.6410924 181.056830
[208,] 105.335738 98.5524738 179.280944
[209,] 104.221403 98.3085533 179.213661
[210,] 103.280753 97.8924620 178.068998
[211,] 102.550153 97.8395338 177.181684
[212,] 102.614671 97.6923930 175.273280
[213,] 101.827870 97.6018089 173.724035
[214,] 100.901485 97.2632299 173.401112
[215,] 100.201867 97.0555162 172.076492
[216,] 99.352513 97.0163491 170.921400
[217,] 98.593180 96.7882731 170.159733
[218,] 97.853335 96.6105136 169.242697
[219,] 97.466315 96.6417549 167.919909
[220,] 96.548521 96.1799468 166.396024
[221,] 95.518471 96.1145910 166.423632
[222,] 94.845758 96.0468343 165.311850
[223,] 94.523744 96.1512785 164.596490
[224,] 93.513806 96.0175943 164.484434
[225,] 92.779019 95.7512629 163.371182
[226,] 92.540550 95.6970929 161.858189
[227,] 92.218382 95.5928066 160.886478
[228,] 91.267649 95.5047679 160.467959
[229,] 90.775148 95.4820795 159.705582
[230,] 90.565735 95.6951957 159.338035
[231,] 89.846931 95.5787396 158.736440
[232,] 89.293953 95.6129115 158.274820
[233,] 88.843813 95.4286769 157.656622
[234,] 88.234318 95.3500693 156.954268
[235,] 87.886669 95.5311066 156.386035
[236,] 87.388184 95.7231352 155.875580
[237,] 87.348655 95.7825825 154.837452
[238,] 86.506522 95.8437068 154.865861
[239,] 86.361072 95.7094117 153.811236
[240,] 85.519379 95.7326171 154.009034
[241,] 85.313088 95.6683310 153.041002
[242,] 84.719943 95.6682763 152.819855
[243,] 84.055191 95.4482427 151.989442
[244,] 83.845519 95.5570044 151.174997
[245,] 83.041307 95.6727176 150.734613
[246,] 82.763247 95.4405915 150.048871
[247,] 82.206473 95.0524500 149.271784
[248,] 81.778340 95.4621719 149.381010
[249,] 81.335401 95.3453500 148.520721
[250,] 81.403587 95.4250837 147.542226
[251,] 80.722578 95.3854522 146.858745
[252,] 80.833523 95.2775211 146.182450
[253,] 80.253653 95.3199663 145.829431
[254,] 79.740133 94.9237396 144.578935
[255,] 79.435585 95.0532062 144.103742
[256,] 78.829549 94.8869302 143.252742
[257,] 78.451354 94.6783548 142.614271
[258,] 77.761259 94.7945615 142.764083
[259,] 77.258806 94.5848047 142.244019
[260,] 77.100983 94.5875330 141.192032
[261,] 76.968030 94.6263475 140.247166
[262,] 76.521581 94.3024009 139.767875
[263,] 75.999684 94.5562247 139.506416
[264,] 75.891508 94.4147972 138.304081
[265,] 75.655621 94.6528889 138.097996
[266,] 75.594235 94.5072096 137.281704
[267,] 75.165432 94.4074324 136.909102
[268,] 74.970210 94.6429197 136.874606
[269,] 74.598544 94.4201499 136.571332
[270,] 74.419679 94.3898679 135.623790
[271,] 74.084147 94.5705727 135.430380
[272,] 73.647409 94.4719941 134.780073
[273,] 73.623847 94.4853537 134.274609
[274,] 73.475858 94.6956620 133.812457
[275,] 73.263236 94.7332009 133.537729
[276,] 73.204971 94.7638937 133.373090
[277,] 72.974472 94.8572602 133.110048
[278,] 72.998243 95.2735880 132.720986
[279,] 72.738265 95.2116584 132.906944
[280,] 72.744121 95.5105281 132.348176
[281,] 72.638573 95.6663775 131.808173
[282,] 72.823548 95.9463536 131.570332
[283,] 72.556854 95.9556729 131.147319
[284,] 72.612351 96.0811926 131.042762
[285,] 72.440545 96.1602553 130.613987
[286,] 72.279979 96.1612475 130.464944
[287,] 72.165040 96.4216479 130.226819
[288,] 71.776503 96.2299095 130.075279
[289,] 71.953473 96.5118170 129.537512
[290,] 71.764427 96.6551087 129.620354
[291,] 71.760978 96.8662763 129.308276
[292,] 71.798116 97.0492060 128.704293
[293,] 71.543434 96.8312613 128.046028
[294,] 71.346362 96.8805494 127.690261
[295,] 71.412662 96.9170203 127.201880
[296,] 71.290508 96.8338554 126.573302
[297,] 71.102466 96.7523767 126.265234
[298,] 70.805117 96.7847298 125.635177
[299,] 70.791130 96.7904069 125.444533
[300,] 70.868267 96.4732788 125.421771
[301,] 70.493065 96.3087908 125.112982
[302,] 70.404410 96.4964430 124.522419
[303,] 70.531671 96.7937439 124.128940
[304,] 70.602699 96.7109476 123.821249
[305,] 70.334441 96.8692391 123.550471
[306,] 70.170046 96.8043253 123.052652
[307,] 70.121549 96.7478605 122.755209
[308,] 70.354025 96.7601385 122.234825
[309,] 70.253305 97.0705875 122.160561
[310,] 70.135896 97.1464845 121.913085
[311,] 69.874392 97.0163407 121.258631
[312,] 69.785899 97.0696455 121.110859
[313,] 70.003710 97.1712884 120.755678
[314,] 69.848822 97.1963879 120.401093
[315,] 69.922598 97.4275375 120.170650
[316,] 69.927257 97.6765706 119.890643
[317,] 69.775766 97.8146695 119.601217
[318,] 69.865381 97.9265590 119.705000
[319,] 69.939432 98.1633412 119.520068
[320,] 70.013134 97.7196310 119.150946
[321,] 70.240042 97.8197485 118.964301
[322,] 70.088034 98.3800204 118.778117
[323,] 70.176711 97.9893529 118.486551
[324,] 70.607859 98.5796317 118.670520
[325,] 70.533625 98.5022533 118.355086
[326,] 70.545206 98.6435829 118.376397
[327,] 70.591806 98.6490788 118.082394
[328,] 70.668438 98.7125502 117.766900
[329,] 70.837114 98.8184388 117.385214
[330,] 70.690716 98.7074804 116.833833
[331,] 70.646488 98.8648709 117.119237
[332,] 70.590418 98.6900362 116.857496
[333,] 70.766485 98.8065552 116.848736
[334,] 70.805570 98.6694389 116.283979
[335,] 70.778672 98.7607651 116.180466
[336,] 70.976094 99.1230099 116.245108
[337,] 70.891753 98.9506982 115.762496
[338,] 70.902807 98.8635402 115.108349
[339,] 70.698457 98.6992138 115.014599
[340,] 70.846168 98.3744105 114.664121
[341,] 70.462085 98.2095985 113.933103
[342,] 70.283082 98.1509037 113.453540
[343,] 70.454126 98.2933285 112.691644
[344,] 70.243580 97.7580172 112.129044
[345,] 70.193647 97.6211087 111.958556
[346,] 70.078341 97.7549259 111.852585
[347,] 70.134841 97.3919161 111.198859
[348,] 69.591741 96.9765731 110.621779
[349,] 69.555315 96.6243819 109.772521
[350,] 69.150043 96.3762015 109.520727
[351,] 69.487061 96.2033717 108.788584
[352,] 68.794998 95.9360967 108.599812
[353,] 68.785639 95.9470727 108.250111
[354,] 68.438332 95.4908172 107.577582
[355,] 68.211680 95.2043990 107.127901
[356,] 67.812791 94.5537628 106.210793
[357,] 67.465485 94.1198868 105.561330
[358,] 67.363219 93.8823817 104.877522
[359,] 67.102348 93.1326387 104.287023
[360,] 66.958672 93.0052612 103.974080
[361,] 67.141048 92.8064252 103.062775
[362,] 66.918342 91.9980705 102.229498
[363,] 66.730119 91.9818140 101.532140
[364,] 66.463332 91.7456215 101.308491
[365,] 66.214728 91.0710338 100.535247
[366,] 65.912155 90.8775234 100.436073
[367,] 65.937461 90.5597000 99.340639
[368,] 65.835298 90.5662078 99.374654
[369,] 65.340586 90.1352563 98.787900
[370,] 65.116993 89.5522571 98.386044
[371,] 65.124216 89.4822619 98.030455
[372,] 64.560702 89.1954701 97.614758
[373,] 64.867071 88.6806681 96.962101
[374,] 64.483320 88.6105928 96.397552
[375,] 64.353096 88.5642197 96.158607
[376,] 64.214050 88.0972009 95.195176
[377,] 63.677964 87.4613327 94.554386
[378,] 63.304845 87.3430649 94.433003
[379,] 63.185703 87.1278510 94.155431
[380,] 63.072980 86.7685313 93.813315
[381,] 62.940424 86.5431471 93.337390
[382,] 62.613294 86.1487267 92.886841
[383,] 62.398085 85.7995405 92.325649
[384,] 62.324354 85.5257102 91.946834
[385,] 62.420648 85.2442565 91.204008
[386,] 62.054359 84.8125722 90.599992
[387,] 61.813623 84.4087462 90.217501
[388,] 61.136849 83.8625053 89.711102
[389,] 61.386586 83.4209744 88.998335
[390,] 60.924325 83.5447727 88.745019
[391,] 60.812772 82.7265793 88.161667
[392,] 60.290445 82.6078116 87.610872
[393,] 60.172295 81.8836750 86.656101
[394,] 59.550861 81.3429669 86.072844
[395,] 59.377899 81.0137760 85.530591
[396,] 59.307290 80.3465812 84.771780
[397,] 58.971026 80.2087495 84.550330
[398,] 58.992973 79.8403035 84.097634
[399,] 58.208248 78.9046891 83.049053
[400,] 58.111814 78.4739719 82.556720
[401,] 57.539510 78.1671024 82.082089
[402,] 57.225432 77.8345172 81.453342
[403,] 56.829731 77.2607806 80.811694
[404,] 56.316228 76.5923070 80.174146
[405,] 56.084056 75.9290886 79.573222
[406,] 55.500846 75.5724531 79.236260
[407,] 55.223778 74.9373075 78.617765
[408,] 54.947377 74.4694775 77.999512
[409,] 54.723625 73.9806539 77.361692
[410,] 54.336448 73.3371910 76.697368
[411,] 53.819250 72.9158745 75.905971
[412,] 53.477569 72.1481178 75.018336
[413,] 53.089046 71.8090646 74.676531
[414,] 52.787577 71.4121376 73.999247
[415,] 52.621226 71.0932569 73.513707
[416,] 51.996956 70.2681192 72.698755
[417,] 51.702334 69.7560198 72.171777
[418,] 51.263349 69.3132509 71.543478
[419,] 50.999568 68.6450940 70.531429
[420,] 50.647057 68.0530806 69.839189
[421,] 50.273899 67.6058682 69.483991
[422,] 49.742837 67.0734402 69.074112
[423,] 49.315085 66.5613991 68.593792
[424,] 48.906282 65.9481205 67.931353
[425,] 48.869655 65.4877590 67.318062
[426,] 48.408860 65.1133096 66.848963
[427,] 48.093356 64.4522535 66.122402
[428,] 47.885292 64.2250814 65.837219
[429,] 47.552243 63.6491926 65.430450
[430,] 47.198983 63.3085580 64.929072
[431,] 46.976763 62.8463555 63.866345
[432,] 46.506963 62.2155484 63.491783
[433,] 46.160509 61.6469987 63.099265
[434,] 45.625284 61.0133060 62.236673
[435,] 45.528323 60.7179898 62.167330
[436,] 45.075531 59.8908052 61.425378
[437,] 44.859211 59.6915714 61.103446
[438,] 44.522991 59.0415037 60.295042
[439,] 44.236413 58.8182428 59.970745
[440,] 43.902635 58.2198016 58.998117
[441,] 43.482127 57.8072981 58.759916
[442,] 43.008733 57.4327706 58.530083
[443,] 42.916562 57.0910084 58.065906
[444,] 42.736564 56.7135270 57.609575
[445,] 42.420203 56.3146550 56.929399
[446,] 42.197105 55.8569975 56.610375
[447,] 41.880534 55.1243086 55.915470
[448,] 41.439165 54.5437864 55.290622
[449,] 41.281269 54.5114432 55.124416
[450,] 40.848898 53.7409114 54.453948
[451,] 40.407669 53.2586410 54.092018
[452,] 40.260678 52.8520817 53.456213
[453,] 39.745399 52.4605903 52.927159
[454,] 39.437811 52.0027810 52.485045
[455,] 39.063678 51.4410683 52.116032
[456,] 38.845305 50.9189721 51.549922
[457,] 38.619696 50.6058922 51.163774
[458,] 38.408089 50.3478571 50.595062
[459,] 38.138368 49.7591114 50.177244
[460,] 37.636884 49.2559570 49.615593
[461,] 37.146041 48.7490251 49.307074
[462,] 36.843877 48.2950388 48.620835
[463,] 36.484139 47.9623060 47.992868
[464,] 36.132799 47.2826583 47.597230
[465,] 35.822218 46.7906017 47.039541
[466,] 35.424274 46.2019437 46.449483
[467,] 35.366665 45.7672959 45.914468
[468,] 34.944337 45.2994888 45.488368
[469,] 34.557196 44.8253570 45.118191
[470,] 34.181615 44.3431641 44.538207
[471,] 33.997883 43.9219635 44.130198
[472,] 33.527514 43.4300813 43.642760
[473,] 33.037279 42.5837760 42.785981
[474,] 32.857912 42.5492558 42.624348
[475,] 32.338975 41.6157093 41.989899
[476,] 32.298649 41.2597636 41.527587
[477,] 31.915716 40.9048877 41.068086
[478,] 31.649197 40.4520438 40.549889
[479,] 31.346764 39.7309124 39.956620
[480,] 30.994531 39.3971521 39.566051
[481,] 30.995270 38.9641738 39.140716
[482,] 30.650329 38.5899307 38.780284
[483,] 30.392964 38.0516385 38.200803
[484,] 30.368529 37.6689670 37.796375
[485,] 30.150901 37.1393331 37.181001
[486,] 29.850110 36.6162946 36.739157
[487,] 29.673507 36.1295034 36.352445
[488,] 29.512182 35.8957968 35.946818
[489,] 29.377477 35.3518461 35.499205
[490,] 29.395687 34.9292181 35.122826
[491,] 29.091301 34.4748099 34.654455
[492,] 29.007689 34.1722687 34.300845
[493,] 28.770854 33.8364482 33.907873
[494,] 28.464079 33.4699600 33.417387
[495,] 28.250202 33.0658648 33.047173
[496,] 28.031538 32.7155852 32.668705
[497,] 27.619892 32.2154810 32.166836
[498,] 27.507217 31.8911028 31.830635
[499,] 27.302721 31.6402086 31.579664
[500,] 27.004787 31.0184236 30.961136
[501,] 26.802900 30.8283357 30.754387
[502,] 26.426886 30.3641909 30.414630
[503,] 26.169907 30.2318988 30.220725
[504,] 25.826648 29.8890106 29.873011
[505,] 25.469211 29.6519594 29.595173
[506,] 25.298590 29.3125703 29.298179
[507,] 24.873730 28.9727597 28.912925
[508,] 24.416003 28.7932364 28.655387
[509,] 24.140214 28.5998858 28.326844
[510,] 23.799464 28.3522355 28.048803
[511,] 23.373861 28.1424942 27.794277
[512,] 22.997167 27.8283860 27.437554
[513,] 22.611955 27.5448211 27.159971
[514,] 22.342379 27.2180932 26.847463
[515,] 22.086065 27.0243366 26.601820
[516,] 21.726938 26.6793888 26.238083
[517,] 21.389398 26.4368822 25.980092
[518,] 21.145042 26.1512298 25.699041
[519,] 20.844910 25.9067654 25.464228
[520,] 20.490598 25.7985359 25.190087
[521,] 20.211713 25.5626125 24.970792
[522,] 20.005238 25.2971877 24.696478
[523,] 19.762824 25.1726801 24.510781
[524,] 19.499280 24.8999935 24.273504
[525,] 19.298056 24.7587017 24.039962
[526,] 19.034981 24.4386384 23.651369
[527,] 18.848639 24.2259704 23.433126
[528,] 18.721624 24.0030799 23.219574
[529,] 18.444675 23.7220301 22.903975
[530,] 18.255543 23.4595183 22.599627
[531,] 18.034231 23.1664772 22.302126
[532,] 17.979610 22.8303775 22.102634
[533,] 17.862578 22.5210602 21.923181
[534,] 17.617786 22.2619123 21.621950
[535,] 17.505179 22.0778776 21.416549
[536,] 17.349889 21.8512729 21.235292
[537,] 17.235709 21.5835744 21.014166
[538,] 17.122858 21.4647233 20.813460
[539,] 16.955003 21.1913191 20.492590
[540,] 16.824548 20.9156469 20.283671
[541,] 16.696310 20.6032265 19.993806
[542,] 16.599919 20.3462238 19.774708
[543,] 16.510182 20.0923079 19.567411
[544,] 16.354441 19.7008154 19.317520
[545,] 16.259869 19.4876605 19.138188
[546,] 16.228117 19.2508215 19.018714
[547,] 16.133068 19.0151594 18.830690
[548,] 16.036073 18.7421263 18.643073
[549,] 15.946948 18.5261823 18.432642
[550,] 15.877955 18.2749890 18.213944
[551,] 15.803750 17.9681025 17.979374
[552,] 15.730516 17.7996062 17.792235
[553,] 15.722362 17.5850027 17.628365
[554,] 15.634058 17.3849822 17.426879
[555,] 15.642303 17.1934073 17.309121
[556,] 15.594689 16.8640576 17.118360
[557,] 15.520770 16.7454197 16.960399
[558,] 15.477818 16.4796340 16.802293
[559,] 15.455418 16.3392477 16.733105
[560,] 15.433770 16.1771071 16.604196
[561,] 15.388948 15.9217118 16.458106
[562,] 15.419322 15.7590117 16.375666
[563,] 15.369182 15.5948802 16.289092
[564,] 15.376938 15.4044170 16.209636
[565,] 15.363430 15.2604665 16.139103
[566,] 15.360895 15.1452977 16.133887
[567,] 15.321499 14.9911914 16.140242
[568,] 15.354468 14.8634264 16.205643
[569,] 15.352377 14.7269768 16.260211
[570,] 15.422172 14.6197279 16.368587
[571,] 15.459469 14.4675093 16.401015
attr(,"class")
[1] "mrfDepth" "mfdmedian"
Cargando paquete requerido: ggplot2
%R -o median_mf
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
def plot_lambdas(lambdas_3, median_mf, n_samples=None):
"""
Plot intensity data for lambdas_3 with median MFD curves
Parameters:
- lambdas_3: numpy array with shape (n_samples, points, wavelengths)
- median_mf: OrderedDict with median MFD curves
- n_samples: number of samples to plot (default: all samples)
"""
if n_samples is None:
n_samples = lambdas_3.shape[0]
# Limit to specified number of samples
lambdas_subset = lambdas_3[:n_samples]
# Create color gray palette
colors = plt.cm.gray(np.linspace(0, 1, n_samples))
plt.figure(figsize=(15, 6))
# Plot each wavelength for samples
for wavelength in range(lambdas_subset.shape[2]):
for sample in range(n_samples):
plt.plot(lambdas_subset[sample, :, wavelength],
color=colors[sample],
linewidth=.5,
alpha=0.3)
# Color palette for median curves
labels = [r"Mediana de $\lambda = 5$" , r"Mediana de $\lambda = 6$", r"Mediana de $\lambda = 7$"]
# Plot median MFD curves
key, curve = next(iter(median_mf.items()))
plt.plot(curve, linewidth=3, label= labels)
plt.title(f'Lambdas Data - First {n_samples} Samples', fontsize=16)
plt.xlabel('Measurement Points')
plt.ylabel('Intensity')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()
# Ejemplo de uso (descomentar y ajustar según tus datos)
# plot_lambdas(lambdas_3, median_mf=median_mf)
# Ejemplo de uso
plot_lambdas(lambdas_3, median_mf=median_mf)
Outliers del proceso funcional multivariado¶
Teniendo un vector de funciones $\mathbf{X}$
el primer paso es la profundidad como, tomando pointwise
$$\text{MFMD} = \int_T D(\mathbf{X} (t), F_{\mathbf{X} (t)})\omega(t) dt$$
Donde $D$ es una profundidad multivariada y MFMD es una profundidad funcional multivariada.
$$ D: \mathcal{R}^p \to \mathcal{R}$$ $$\text{MFMD}: \mathcal{L}^{2p}_T \to \mathcal{R}_{[0, 1]}$$
Se define una función de outlyingness como $$ o: \mathcal{R}^p \to [0 , \infty)$$
La profundidad funcional multivariada se puede definir con base a la función de outlyingness (atipicidad) como
$$ D(\mathbf{X} (t), F_{\mathbf{X} (t)}) = [1 + o(\mathbf{X} (t), F_{\mathbf{X} (t)})]^{-1}$$ Despejando:
$$ o(\mathbf{X} (t), F_{\mathbf{X} (t)}) = [D(\mathbf{X} (t), F_{\mathbf{X} (t)})^{-1} - 1] $$
Ahora se define la función O por medio de la O:
$$ O((\mathbf{X} (t), F_{\mathbf{X} (t)})) = o((\mathbf{X} (t), F_{\mathbf{X} (t)}))\vec{v}((\mathbf{X} (t), F_{\mathbf{X} (t)}))$$
$$ O: \mathcal{R}^p \to \mathcal{R}^p$$
Donde $\vec{v}(\mathbf{y})$ es un vector unitario en la dirección de $\mathbf{y}$
$$ \vec{v}(\mathbf{y}) = \frac{\mathbf{y} - \tilde{\mathbf{y}}}{||\mathbf{y} - \tilde{\mathbf{y}}||} = \text{Signo del vector y} = Sign(\mathbf{y})$$
Anomalía Direccional Funcional¶
Se define un proceso estocástico $ X : I \to \mathbb{R}^p $ que toma valores en un espacio de funciones continuas $ C(I, \mathbb{R}^p) $, con distribución de probabilidad $ F_X $. Para medir la anomalía direccional de $ X $, se utilizan las siguientes métricas:
Anomalía direccional funcional ($ FO $): Representa la anomalía total de $ X $, acumulando la magnitud de la anomalía direccional a lo largo del intervalo $ I $: $$ FO(X, F_X) = \int_I \|O(X(t), F_{X(t)})\|^2 w(t) \, dt $$ Aquí, $ O(X(t), F_{X(t)}) $ es la anomalía direccional en el tiempo $ t $, y $ w(t) $ es una función de peso.
Anomalía direccional media ($ MO $): Evalúa la posición promedio de $ X $ en relación con la curva central, considerando tanto la distancia como la dirección: $$ MO(X, F_X) = \int_I O(X(t), F_{X(t)}) w(t) \, dt $$ $$ =\int_I\left[ f , f , \dots, f \right]$$ $$ =\left[ \int_I , \int_I , \dots, \int_I \right] \in \mathcal{R}^p$$
- Note que MO devuelve un vector, (se esta integrando un vector), las funciones f pertenecen a $\mathcal{R}^p$
- La norma de MO da la la magnitud de la atipicidad que son outliers de magnitud.
- Variación de la anomalía direccional ($ VO $): Mide las variaciones de $ O(X(t), F_{X(t)}) $ en términos de norma y dirección a lo largo del intervalo: $$ VO(X, F_X) = \int_I \|O(X(t), F_{X(t)}) - MO(X, F_X)\|^2 w(t) \, dt $$
- Note que VO es un número.
- Determina los outliners de forma.
Se puede definir la distancia de mahalanobis robusta como
$$ \vec{y} = (\underbrace{\vec{MO}}_{\mathcal{R}^p}, \underbrace{VO}_{\mathcal{R}}) \in \mathcal{R}^{p +1}$$
La función de covarianza cambia como se usa MCD (minimum covariance determinant) para encontrar la covarianza robusta.
$$ \frac{1}{h} \sum_{i=1}^h (\vec{y}_i - \bar{\vec{y}})(\vec{y}_i - \bar{\vec{y}})^T$$
Distribución y Detección de Valores Atípicos en Curvas Generadas por Procesos Gaussianos¶
A través de estudios numéricos, se encontró que la distribución de $ Y_{k,n} = (MO^T_{T_k,n}, VO_{T_k,n})^T $ se puede aproximar bien mediante una distribución normal de dimensión $ p+1 $ cuando $ X $ es generado por un proceso estacionario gaussiano de $ p $ dimensiones.
Detección de Valores Atípicos:
- Se utilizó la distancia de Mahalanobis robusta ($ RMD^2 $) junto con los estimadores mínimos del determinante de covarianza (MCD)
Procedimiento para Detectar Valores Atípicos¶
Cálculo de la distancia de Mahalanobis robusta: $$ RMD^2 (Y_{k,n}, \bar{Y}^*_{k,n,J}) = (Y_{k,n} - \bar{Y}^*_{k,n,J})^T S^*_{k,n,J}{}^{-1}(Y_{k,n} - \bar{Y}^*_{k,n,J}), $$ donde:
- $ J $: grupo de $ h $ puntos que minimizan el determinante de la matriz de covarianza.
- $ \bar{Y}^*_{k,n,J} = h^{-1} \sum_{i \in J} Y_{k,n,i} $.
- $ S^*_{k,n,J} = h^{-1} \sum_{i \in J} (Y_{k,n,i} - \bar{Y}^*_{k,n,J})(Y_{k,n,i} - \bar{Y}^*_{k,n,J})^T $.
- $ h $: submuestra que controla la robustez del método.
Aproximación de la cola de la distribución: Según Hardin y Rocke (2005), la cola de $ RMD^2 $ sigue: $$ \frac{c(m - p)}{m(p + 1)} RMD^2 (Y_{k,n}, \bar{Y}^*_{k,n,J}) \sim F_{p+1,m-p}, $$ donde:
- $ c $ y $ m $: parámetros para la distribución $ F $ y el factor de escala, calculados mediante simulación.
- El valor de corte ($ C $) es el cuantil $ \alpha $ ($ \alpha = 0.993 $) de $ F_{p+1,m-p} $.
Identificación de valores atípicos: Una curva se clasifica como atípica si: $$ \frac{c(m - p)}{m(p + 1)} RMD^2 (Y_{k,n}, \bar{Y}^*_{k,n,J}) > C. $$
Usos Adicionales del $ RMD $¶
- Medida de centralidad para definir:
- La mediana y la región central de las curvas.
- La función de media recortada.
- El boxplot funcional.
lambdas_3.shape
(268, 571, 3)
%R -i lambdas_3
%%R
library("fdaoutlier")
# Show mr the dimensions of lambdas_3
print(dim(lambdas_3))
directional_outliners <- dir_out(
lambdas_3,
data_depth = "random_projections",
n_projections = 200L,
seed = NULL,
return_distance = TRUE,
return_dir_matrix = TRUE
)
[1] 268 571 3
Con las distancias es necesario hacer re muestreo de n, con 1000 veces ajusto una distribución y los que esten fuera del quantil 99.3
%R -o directional_outliners
# show me the components of the directional_outliners
display(directional_outliners.keys())
display(directional_outliners['mean_outlyingness'].shape)
display(directional_outliners['var_outlyingness'].shape)
display(directional_outliners['ms_matrix'].shape)
display(directional_outliners['distance'].shape)
(np.str_('mean_outlyingness'),
np.str_('var_outlyingness'),
np.str_('distance'),
np.str_('ms_matrix'),
np.str_('mcd_obj'),
np.str_('dirout_matrix'))
(268, 3)
(268,)
(268, 4)
(268,)
distance = directional_outliners['distance']
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
def bootstrap_percentile_distribution(data, n_resamples=1000, percentile=99.3):
"""
Realiza remuestreo bootstrap, grafica la distribución empírica
y marca el percentil especificado.
Parámetros:
-----------
data : list or numpy.ndarray
Lista o array de datos originales
n_resamples : int, opcional (por defecto 1000)
Número de remuestreos bootstrap
percentile : float, opcional (por defecto 99.7)
Percentil para determinar valores extremos
Retorna:
--------
tuple
- Array de índices de valores originales que superan el percentil
- Percentil calculado
"""
# Convertir los datos a un array de NumPy
data = np.array(data)
# Número de muestras
n = len(data)
# Almacenar resultados de los remuestreos
bootstrap_samples = np.zeros((n_resamples, n))
# Realizar remuestreo bootstrap
for i in range(n_resamples):
# Remuestreo con reemplazo
bootstrap_samples[i] = np.random.choice(data, size=n, replace=True)
# Calcular el percentil para cada remuestreo
bootstrap_percentiles = np.percentile(bootstrap_samples, percentile, axis=1)
# Calcular el cuantil de corte
quantil_corte = np.percentile(bootstrap_percentiles, percentile)
# Encontrar los índices de los valores originales por encima del percentil
extreme_indices = np.where(data > quantil_corte)[0]
# Configurar el estilo de la gráfica
plt.figure(figsize=(10, 6))
sns.set_style('whitegrid')
# Graficar la distribución empírica de los remuestreos
sns.histplot(bootstrap_percentiles, kde=True, color='blue', alpha=0.7)
# Añadir línea vertical roja punteada para el cuantil de corte
plt.axvline(x=quantil_corte, color='red', linestyle='--',
label=f'Cuantil {percentile}%: {quantil_corte:.4f}')
plt.title('Distribución Empírica de Percentiles Bootstrap')
plt.xlabel('Distancia robusta de Mahalanobis')
plt.ylabel('Frecuencia')
plt.legend()
plt.tight_layout()
plt.show()
return extreme_indices, quantil_corte
indices_extremos, cuantil = bootstrap_percentile_distribution(distance)
print("Índices de valores extremos:", indices_extremos)
print(f"Cuantil de corte (99.3): {cuantil}")
print("Valores extremos:", distance[indices_extremos])
Índices de valores extremos: [] Cuantil de corte (99.3): 2443.7431231536584 Valores extremos: []
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
def estimate_distribution_parameters(distances,
n_bootstrap: int = 1000,
alpha: float = 0.993,
seed: int = None):
"""
Estimate m and c parameters for the F-distribution based approach
Parameters:
-----------
distances : np.ndarray or list
Input distances to analyze
n_bootstrap : int, optional
Number of bootstrap iterations
alpha : float, optional
Significance level for outlier cutoff
seed : int, optional
Random seed for reproducibility
Returns:
--------
Tuple of (estimated m, estimated c, standard deviation of c)
"""
# Convert to numpy array if not already
distances = np.asarray(distances)
flat_distances = distances.flatten()
if seed is not None:
np.random.seed(seed)
# Number of features (p)
p = 1 # Assuming univariate case
# Prepare arrays to store bootstrap results
m_estimates = np.zeros(n_bootstrap)
c_estimates = np.zeros(n_bootstrap)
# Bootstrapping loop
for i in range(n_bootstrap):
# Resample distances with replacement
bootstrap_sample = np.random.choice(flat_distances, size=len(flat_distances), replace=True)
# Try different m values
m_candidates = np.linspace(len(flat_distances)/2, len(flat_distances)*2, 20)
best_error = np.inf
best_m = len(flat_distances)
best_c = 1.0
for m_candidate in m_candidates:
# Scale distances according to the formula
scaled_distances = bootstrap_sample * (m_candidate - p) / (m_candidate * (p + 1))
try:
# Theoretical F-distribution quantile
theoretical_quantile = stats.f.ppf(alpha, p + 1, m_candidate - p)
# Empirical quantile
empirical_quantile = np.quantile(scaled_distances, alpha)
# Calculate error
error = abs(theoretical_quantile - empirical_quantile)
# Update best parameters if error is smaller
if error < best_error:
best_error = error
best_m = m_candidate
best_c = theoretical_quantile / empirical_quantile
except Exception:
continue
# Store the best estimates
m_estimates[i] = best_m
c_estimates[i] = best_c
# Remove any zero or invalid estimates
m_estimates = m_estimates[m_estimates > 0]
c_estimates = c_estimates[c_estimates > 0]
# Compute statistics
mean_m = np.mean(m_estimates)
mean_c = np.mean(c_estimates)
std_c = np.std(c_estimates)
return mean_m, mean_c, std_c
def detect_outliers(distances,
alpha: float = 0.993,
bootstrap_m_c: bool = True,
n_bootstrap: int = 1000,
plot: bool = True):
"""
Detect outliers using robust distance method
Parameters:
-----------
distances : np.ndarray or list
Input distances to analyze
alpha : float, optional
Significance level for outlier cutoff
bootstrap_m_c : bool, optional
Whether to estimate m and c via bootstrapping
n_bootstrap : int, optional
Number of bootstrap iterations
plot : bool, optional
Whether to generate diagnostic plots
Returns:
--------
List of outlier indices
"""
# Convert to numpy array if not already
distances = np.asarray(distances)
flat_distances = distances.flatten()
# Number of features (p)
p = 1 # Assuming univariate case
# Estimate m and c via bootstrapping if requested
if bootstrap_m_c:
m, c_factor, c_std = estimate_distribution_parameters(
distances,
n_bootstrap=n_bootstrap,
alpha=alpha
)
else:
# Default values if not bootstrapping
m = len(flat_distances)
c_factor = 1.0
c_std = None
# Scale robust distances
scaled_distances = c_factor * (m - p) / (m * (p + 1)) * flat_distances
# Cutoff value from F-distribution
cutoff_value = stats.f.ppf(alpha, p + 1, m - p)
# Detect outliers
is_outlier = scaled_distances > cutoff_value
# Create diagnostic plots if requested
if plot:
plt.figure(figsize=(15, 5))
# Original Distances
plt.subplot(1, 3, 1)
plt.hist(flat_distances, bins=30, density=True, alpha=0.7)
plt.title('Original Distances')
plt.xlabel('Distance')
# Scaled Distances
plt.subplot(1, 3, 2)
plt.hist(scaled_distances, bins=30, density=True, alpha=0.7)
plt.title('Scaled Distances')
plt.axvline(cutoff_value, color='r', linestyle='--', label='Cutoff')
plt.xlabel('Scaled Distance')
plt.legend()
# Box Plot of Distances
plt.subplot(1, 3, 3)
plt.boxplot(flat_distances)
plt.title('Boxplot of Distances')
plt.xlabel('Distance')
plt.tight_layout()
plt.show()
# Q-Q Plot
plt.figure(figsize=(10, 5))
stats.probplot(flat_distances, dist="f", sparams=(p+1, m-p), plot=plt)
plt.title('Q-Q Plot of Distances')
plt.show()
# Print diagnostic information
print(f"Estimated m: {m:.4f}")
print(f"Estimated c factor: {c_factor:.4f}")
if c_std is not None:
print(f"Standard deviation of c: {c_std:.4f}")
print(f"Cutoff value: {cutoff_value:.4f}")
print(f"Number of outliers: {sum(is_outlier)}")
# Return indices of outliers
return np.where(is_outlier)[0].tolist()
print("\n--- Without Bootstrapping M and C ---")
outliers_without_bootstrap = detect_outliers(distance, bootstrap_m_c=False)
print("\n--- With Bootstrapping M and C ---")
outliers_with_bootstrap = detect_outliers(distance, bootstrap_m_c=True)
print(f"Los outliers detectados son: {outliers_with_bootstrap}")
# con retorno
--- Without Bootstrapping M and C ---
Estimated m: 268.0000 Estimated c factor: 1.0000 Cutoff value: 5.0552 Number of outliers: 103 --- With Bootstrapping M and C ---
Estimated m: 134.0000 Estimated c factor: 0.0118 Standard deviation of c: 0.0050 Cutoff value: 5.1516 Number of outliers: 3 Los outliers detectados son: [70, 130, 267]
Tareas Taller 1¶
Profundidad de variación total y similitud de forma¶
import os
import numpy as np
import skfda
from skfda.representation.basis import BSplineBasis
from skfda.preprocessing.smoothing import BasisSmoother
import matplotlib.pyplot as plt
from tqdm import tqdm
def smooth_functional_data(
data,
wavelengths,
dominio,
n_basis=65,
best_param=1e-3,
rango_curvas=[None, None]
):
# Determinar el rango de curvas a procesar
inicio = rango_curvas[0] if rango_curvas[0] is not None else 0
fin = rango_curvas[1] if rango_curvas[1] is not None else data.shape[0]
# Inicializar array para almacenar datos suavizados
smoothed_data = np.zeros_like(data[inicio:fin])
# Crear una barra de progreso
with tqdm(total=(fin-inicio)*len(wavelengths), desc="Suavizando muestras") as pbar:
# Para cada longitud de onda
for wave_idx, _ in wavelengths.items():
# Procesar muestras para esta longitud de onda
for i, muestra in enumerate(range(inicio, fin)):
# Obtener los datos para esta muestra y longitud de onda
X = data[muestra, :, wave_idx].flatten()
# Crear una base B-Spline
basis = BSplineBasis(domain_range=(275, 560), n_basis=n_basis)
# Crear functional data grid
fd = skfda.FDataGrid(data_matrix=X, grid_points=dominio)
# Suavizado penalizado
smoother = BasisSmoother(basis, smoothing_parameter=best_param)
fd_smoothed = smoother.fit_transform(fd)
# Guardar los datos suavizados (flatten to ensure correct shape)
smoothed_data[i, :, wave_idx] = fd_smoothed(dominio)[0].flatten()
# Actualizar la barra de progreso
pbar.update(1)
return smoothed_data
init = 5*268
end = 6*268
# Ejemplo de uso
smoothed_result = smooth_functional_data(data, wavelengths, dominio)
Suavizando muestras: 100%|██████████| 1876/1876 [05:30<00:00, 5.68it/s]
print(class_outliners.keys())
# Obtener los outliers de forma y magnitud
tvd = class_outliners['tvd']
mss = list(class_outliners['mss'])
shp_out = list(class_outliners['shape_outliers'])
# has que shp_out no tenga componentes del tipo np.int32
shp_out = [int(i) for i in shp_out]
mg_out = list(class_outliners['magnitude_outliers'])
(np.str_('outliers'), np.str_('shape_outliers'), np.str_('magnitude_outliers'), np.str_('tvd'), np.str_('mss'))
1. Boxplot de MSS¶
import matplotlib.pyplot as plt
import numpy as np
def boxplot_mss(data, f=1.5):
"""
Crea un boxplot y detecta los outliers de una lista de valores.
Args:
data (list): Lista de valores numéricos
f (float): Factor para calcular los límites de los outliers (default: 1.5)
Returns:
list: Lista con los índices de los outliers encontrados
"""
# Crear la figura
plt.figure(figsize=(8, 6))
# Crear el boxplot
bp = plt.boxplot(data, patch_artist=False)
# Calcular estadísticas
Q1 = np.percentile(data, 25)
Q3 = np.percentile(data, 75)
IQR = Q3 - Q1
lower_bound = Q1 - f * IQR
upper_bound = Q3 + f * IQR
# Identificar índices de outliers
outlier_indices = [i for i, x in enumerate(data) if x < lower_bound or x > upper_bound]
outlier_values = [data[i] for i in outlier_indices]
# Plotear los outliers en color naranja
if len(outlier_values) > 0:
plt.plot([1] * len(outlier_values), outlier_values, 'o', color='orange', label='Outliers')
# Configurar el gráfico
plt.title('Shape outlyingness Boxplot')
plt.grid(True, linestyle='--', alpha=0.7)
# Añadir texto con los límites
plt.figtext(0.02, 0.02,
f'Límite inferior: {lower_bound:.2f}\n'
f'Límite superior: {upper_bound:.2f}',
fontsize=8)
# Mostrar el gráfico
plt.show()
return outlier_indices
shape_out = boxplot_mss(mss, 3)
# a todos los valores de shape_out súmalele 1
shape_out_correction = [i+1 for i in shape_out]
print(f"Índices de los outliers (py): {shape_out_correction}")
print(f"Índices de los outliers (R ): {shp_out}")
Índices de los outliers (py): [104, 109, 125, 129, 138, 140, 143, 145, 147, 149, 150, 151, 152, 153, 154, 155, 159, 160, 162, 163, 164, 167] Índices de los outliers (R ): [104, 109, 125, 129, 138, 140, 143, 145, 147, 149, 150, 151, 152, 153, 154, 155, 159, 160, 162, 163, 164, 167]
2. functional boxplot TVD¶
def functional_boxplot_tvd(X, shape_out, F=1.5):
"""
Create functional boxplot with outliers
"""
# Filtra las filas tal que NO esten en shape_out
X_out = X
fd = skfda.FDataGrid(data_matrix=X_out, grid_points=dominio, dataset_name=f'Functional Boxplot TVD')
fdBoxplot = Boxplot(fd, factor=F)
fdBoxplot.plot()
return fdBoxplot
fboxplot = functional_boxplot_tvd(X, shape_out, 1.5)
outliers = fboxplot.outliers
# get the index where the outliers are
index_out = list(np.where(outliers)[0])
magnitude_out = [int(x) for x in index_out]
index_out = [x + 1 for x in magnitude_out]
mg_out = [int(x) for x in mg_out]
print(f"Índices de los outliers (py): {index_out}")
print(f"Índices de los outliers (R ): {mg_out}")
Índices de los outliers (py): [10, 71, 131] Índices de los outliers (R ): [10, 71, 131]
import numpy as np
import matplotlib.pyplot as plt
def plot_kind_class_outliners(magnitude_outliers, shape_outliers, X):
"""
Plot outliers and non-outliers with different colors and labels.
Shows only 3 shape outliers and all other curves in grey with alpha=0.5
Args:
magnitude_outliers (list): List of indices for magnitude outliers
shape_outliers (list): List of indices for shape outliers
X (numpy.ndarray): Array with all curves data
"""
# Convierte los índices de outliers a conjuntos para una búsqueda eficiente
magnitude_outliers_set = set(magnitude_outliers)
shape_outliers_set = set(shape_outliers)
# Crear máscaras para diferentes tipos de outliers
magnitude_outliers_mask = np.array([i in magnitude_outliers_set for i in range(len(X))])
shape_outliers_mask = np.array([i in shape_outliers_set for i in range(len(X))])
# Encuentra los magnitude outliers en X
magnitude_curves = X[magnitude_outliers_mask, :]
# Encuentra los shape outliers en X (limitado a 3)
shape_curves = X[shape_outliers_mask, :]
shape_curves = shape_curves[:3] if len(shape_curves) > 3 else shape_curves
# Encuentra las curvas que no son outliers
non_outliers_mask = ~(magnitude_outliers_mask | shape_outliers_mask)
non_outliers = X[non_outliers_mask, :]
# Grafica los tipos de curvas
plt.figure(figsize=(12, 6))
# Primero grafica todas las curvas no outliers en gris
for i in range(non_outliers.shape[0]):
plt.plot(non_outliers[i], color='gray', alpha=0.5, label='_')
# Grafica los magnitude outliers
for i in range(magnitude_curves.shape[0]):
plt.plot(magnitude_curves[i], color='red', alpha=0.7,
label=('' if i==0 else '_') + 'Outliers de Magnitud')
# Grafica solo 3 shape outliers
for i in range(len(shape_curves)):
plt.plot(shape_curves[i], color='blue', alpha=0.7,
label=('' if i==0 else '_') + 'Outliers de Forma')
plt.title('Outliers de Magnitud y Forma')
plt.xlabel('Longitud de onda (nm)')
plt.ylabel('Intensidad')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# Imprime estadísticas de outliers
print(f"Outliers de magnitud: {len(magnitude_outliers)}")
print(f"Outliers de forma: {len(shape_outliers)}")
print(f"Shape outliers mostrados: {min(3, len(shape_outliers))}")
# Ejemplo de uso:
plot_kind_class_outliners(magnitude_out, shape_out, smoothed_result[:, :, wavelength_idx])
Outliers de magnitud: 3 Outliers de forma: 22 Shape outliers mostrados: 3
import numpy as np
import matplotlib.pyplot as plt
def shape_outliers(magnitude_outliers, shape_outliers, X):
"""
Plot outliers and non-outliers with different colors and labels.
Shows only 3 shape outliers and all other curves in grey with alpha=0.5
Args:
magnitude_outliers (list): List of indices for magnitude outliers
shape_outliers (list): List of indices for shape outliers
X (numpy.ndarray): Array with all curves data
"""
# Convierte los índices de outliers a conjuntos para una búsqueda eficiente
magnitude_outliers_set = set(magnitude_outliers)
shape_outliers_set = set(shape_outliers)
# Crear máscaras para diferentes tipos de outliers
magnitude_outliers_mask = np.array([i in magnitude_outliers_set for i in range(len(X))])
shape_outliers_mask = np.array([i in shape_outliers_set for i in range(len(X))])
# Encuentra los magnitude outliers en X
magnitude_curves = X[magnitude_outliers_mask, :]
# Encuentra los shape outliers en X (limitado a 3)
shape_curves = X[shape_outliers_mask, :]
shape_curves = shape_curves[19:] if len(shape_curves) > 3 else shape_curves
# Encuentra las curvas que no son outliers
non_outliers_mask = ~(magnitude_outliers_mask | shape_outliers_mask)
non_outliers = X[non_outliers_mask, :]
# Grafica los tipos de curvas
plt.figure(figsize=(12, 6))
# Primero grafica todas las curvas no outliers en gris
for i in range(25):
plt.plot(non_outliers[i], color='gray', alpha=1, label='_')
# Grafica solo 3 shape outliers
for i in range(len(shape_curves)):
plt.plot(shape_curves[i], color='red', alpha=1,
label=('' if i==0 else '_') + 'Outliers de Forma')
plt.title('Outliers de Magnitud y Forma')
plt.xlabel('Longitud de onda (nm)')
plt.ylabel('Intensidad')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# Imprime estadísticas de outliers
print(f"Outliers de magnitud: {len(magnitude_outliers)}")
print(f"Outliers de forma: {len(shape_outliers)}")
print(f"Shape outliers mostrados: {min(3, len(shape_outliers))}")
# Ejemplo de uso:
shape_outliers(magnitude_out, shape_out, smoothed_result[:, :, wavelength_idx])
Outliers de magnitud: 3 Outliers de forma: 22 Shape outliers mostrados: 3